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TECHNICAL  REPORT  SUMMARY 


In  this  contract,  we  set  out  to  find  practical  methods  of  using  the  pulse 
shape  of  observed  seismic  body  waves  to  infer  as  much  as  possible  about  the 
motions  which  must  have  occurred  at  the  seismic  source  itself.  Previous 
attempts  on  this  problem  had  given  insufficient  attention  to  the  constraints 
which  follow  from  requiring  a  physically  plausible  stress  on  the  fault  plane. 

We  believe  we  have  made  substantial  progress  towards  the  solution  of  this  pro¬ 
blem.  Full  details  of  what  we  have  done  are  contained  in  papers  already  pub¬ 
lished  or  now  in  preparation.  It  is  likely  that  a  total  of  six  papers  will 
result  from  the  support  provided  by  this  contract.  They  are  as  follows: 

(i)  Detailed  Spectral  Analysis  of  Two  Small  New  York  Earthquakes , 
by  J.  Boatwright,  published  in  the  Bulletin  of  the  Seismological  Society  of 
America,  68,  1117-1131,  1978.  This  work  was  described  in  full  in  a  previous 
semi-annual  report.  It  gave  a  complete  account  of  the  inversion  of  amplitude 
spectra  of  strong-motion  accelerograms  obtained  at  epicentral  distances  of  only 
about  1  kilometer,  to  obtain  source  dimensions  and  stress  drops  for  two  events 
in  the  eastern  U.S. 

(ii)  A  Spectral  Theory  for  Circular  Seismic  Sources ;  Simple  Esti¬ 
mates  of  Source  Dimension ,  Effective  Stress  and  Radiated  Seisittic  Energy ,  by 
J.  Boatwright,  submitted  for  publication  in  the  Bulletin  of  the  Seismological 
Society  of  America.  This  work  was  described  in  full  in  a  previous  semi-annual 
report.  It  described  the  pulse  shapes  to  be  expected  for  a  fault  rupture  which 
initiated  at  a  point  and  subsequently  grew  as  an  expanding  circle.  A  variety 
of  stopping  mechanisms  were  discussed,  and  it  was  shown  that  a  useful  method 

of  data  processing  involved  working  with  the  square  of  the  observed  particle 
velocity  recorded  at  a  given  station. 

(iii)  Quasi-Dynamic  Models  of  Simple  Earthquakes  and  the  Impli¬ 
cations  of  Energy  Flux  Pulse  Shapes  as  Modelling  Constraints,  by  J.  Boatwright, 
now  in  final  stages  of  preparation  for  submission  probably  to  the  Bulletin  of 
the  Seismological  Society  of  America.  A  copy  of  this  manuscript  is  included 

in  this  report.  A  model  of  fault  slip  that  is  satisfactory  during  stages  of 
rupture,  healing  and  stopping  is  discussed  in  terms  of  the  far -field  pulse 
shapes  it  will  generate.  An  application  is  given  for  two  earthquakes  which 
occurred  in  Alaska.  Estimates  are  given  of  rupture  velocity  and  stress  re¬ 
lease. 


Civ)  Investigations  of  Two  High  Stress-Drop  Earthquakes  in  the 
Shumagin  Seismic  Gap,  Alaska ,  by  L.  House  and  J.  Boatwright.  It  is  currently 
intended  that  this  paper  will  be  submitted  for  publication  in  the  Journal  of 
Geophysical  Research.  A  copy  of  the  manuscript  is  given  in  this  report.  It 
is  inferred  that  the  two  earthquakes  analysed  had  stress-drops  in  excess  of 
500  bars. 


Cv)  Body  Wave  Analysis  of  the  St.  Elias  Earthquake,  section  written 
by  J.  Boatwright  for  inclusion  (with  contributions  from  other  authors)  in  a 
major  paper  on  this  large  and  recent  event.  From  the  far -field  pulse  shapes, 
it  was  inferred  that  three  separate  sub-events  could  be  distinguished  within 
this  rupture.  A  copy  of  the  manuscript  In  Included  in  this  report. 
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(vi)  Elementary  Solutions  to  Lamb's  Problem  for  a  Point  Source  and 
Their  Relevance  to  the  Study  of  Spontaneous  Crack  Propagation  in  Three  Dimen¬ 
sions, ,  by  Paul  G.  Richards,  accepted  for  publication  in  the  August  1979  issue 
of  the  Bulletin  of  the  Seismological  Society  of  America.  In  order  to  determine 
the  slip  occurring  on  a  rupturing  fault  plane,  the  dynamic  consideration  of 
effects  of  initial  stress  and  varying  fault  strength  can  be  handled  if  a 
certain  fundamental  problem  in  elasticity  theory  is  solved  in  a  fashion  that 
admits  rapid  computation  of  the  solution.  This  paper  (a  copy  is  included  in 
this  report)  gives  such  a  solution,  and  discusses  its  significance  in  the 
analysis  of  rupturing  fault  planes. 

The  four  sections  that  follow  are  the  expanded  accounts  of  (Hi)  ,  (iv) , 

(v)  and  (vi)  above. 


QUAS I -DYNAMIC  MODELS  OF  SIMPLE  EARTHQUAKES  AMD  THE 
IMPLICATIONS  OF  ENERGY  FLUX  PULSE  SHAPES 
AS  MODELLING  CONSTRAINTS1 

John  Boatwright 

Lamont-Doherty  Geological  Observatory  and 
Department  of  Geological  Sciences 
of  Columbia  University 
Palisades,  New  York  10964 


ABSTRACT 


We  have  designed  a  code  by  which  one  can  compute  the  far-fleld 
body  wave  displacement  and  energy  flux  pulse  shapes  from  a  series  of 
"quasi-dynamic"  models  of  rupture.  The  Integration  over  the  (kine¬ 
matic)  slip  velocity  is  calculated  on  a  radial  grid  using  the 
Fraunhofer  approximation.  The  specification  of  the  slip  velocity  on 
the  grid  is^ derived  from  theoretical  and  finite  difference  solutions 
for  the  mixed  boundary  value  problem  of  a  3-D  frictional  rupture 
model.  The  analytic  form  for  the  slip  velocity  is  naturally  divided 
into  two  phases:  the  rupture  growth,  in  which  the  slip  distribution 
Is  self -similar  and  elliptical,  and  the  healing,  during  which  the  slip 
velocity,  multiplied  by  a  linearly  decreasing  factor,  goes  to  zero. 

The  arrival  of  a  P-wave  stopping  phase,  generated  by  the  stopping  events 
which  determine  the  fault  perimeter,  determines  the  onset  of  the  healing 
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The  form  of  the  healing  factor,  applied  to  the  rapture  front,  produces 
a  realistic  stopping  of  the  rupture,  enabling  one  to  model  directly 
a  vide  range  of  high-frequency  body  wave  spectra. 

These  quasi-dynamic  models  yield  radiation  efficiencies  consistent 
with  the  theoretical  results  of  Kostrov  0-974} ,  so  that  the  computed 
energy  flux  pulse  shapes  and  the  time-integrated  energy  flux  may  be 
used  as  constraints  in  the  modelling  of  simple  ruptures.  In  particular, 
the  distinct  variation  of  the  energy  flux  pulse  shapes  provides  seismo¬ 
logists  with  a  useful  model  discriminant,  with  Implications  for  the 
determination  of  both  rupture  grovth  and  stopping  behavior  for  multiply 
recorded  earthquakes. 

We  have  applied  this  waveform  modelling  approach  to  the  analysis 
of  two  high  stress  drop  earthquakes  which  occurred  In  the  Shumagln 
Islands,  Alaska,  on  May  6,  1974.  We  analysed  both  short  period  WWSSN 
data  and  strong  motion  accelerograph  data  obtained  from  an  SMA1  at 
Sand  Point.  The  energy  flux  modelling  is  shown  to  provide  an  estimate 
of  the  rupture  velocity  and  thereby  establish  closer  bounds  for  the 
estimates  of  stress  release. 
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INTRODUCTION 


The  forward  problem  of  modelling  seismic  sources  by  matching  the 
waveforms  observed  at  a  restricted  set  of  stations  Is  strongly  non¬ 
unique,  as  has  been  shown  by  a  number  of  authors,  notably  Anderson  and 
Richards  0.976}  and  Boore  and  Stlerman  (1976).  These  and  other  papers 
have  established  the  possibility  of  fitting  the  observed  waveforms 
equally  well  with  different  source  models  incorporating  a  wide  range  of 
prescribed  slip-functions.  Analysis  of  the  inverse  problem  of  seismic 
source  theory  (Kostrov,  1975}  indicates  that  this  non-uniqueness  results 
from  the  Inherent  nature  of  the  seismic  observations,  in  particular  the 
integral  nature  of  the  source  representation,  rather  than  from  inade¬ 
quate  station  coverage. 

While  the  appropriate  slip  function  cannot  be  resolved  from  seismic 
observations,  the  final  source  models  determined  using  different  slip 
functions  can  differ  substantially,  particularly  in  their  estimates  of 
rupture  velocity  and  source  dimension.  Since  both  the  usual  modelling 
estimates  of  stress  released  in  tectonic  earthquakes  (l.e.,  the  dynamic 
stress  drop,  t  ,  and  the  static  stress  drop,  A a)  depend  strongly  on 
these  parameters  (Boatwright,  1979),  the  problem  of  choosing  an  appro¬ 
priate  slip  function  is  of  obvious  seismologlcal  significance.  It  is 
this  choice  which  this  paper  seeks  to  restrict  in  a  heuristic  fashion, 
through  the  use  of  dynamically  consistent  rupture  models. 

In  order  to  determine  a  suitable  source  model,  or  claas  of  source 
models,  It  is  necessary  to  consider  both  a  general  fracture  model  and 
an  appropriate  set  of  initial  conditions.  Our  fracture  model  Is 


derived  from  a  frictional  theory  of  rupture,  l.e.,  the  rupture  Is 
modelled  as  a  growing  region  of  stress  relaxation,  where  the  stress 
acting  on  the  (unhealed)  rupture  area  Is  specified  to  be  the  dynamic 
frictional  stress.  We  presume  the  Initial  loading  stress  to  be  approx* 
lmately  constant  or  to  decrease  away  from  the  rupture  origin,  and  the 
final  rupture  area  to  be  simple  and  planar .  These  assumptions  restrict 
the  application  of  our  rupture  model  to  small  and  moderate  sized  earth¬ 
quakes,  as  large  earthquakes  may  have  strongly  complex  ruptures 
GCanamorl  and  Stewart,  1977) . 

To  complete  our  model  description,  we  need  to  assume  a  fracture 
criterion  and  a  suitable  distribution  of  fracture  strength,  which  will 
then  determine  the  motion  of  the  crack  tip.  For  simplicity,  however, 
the  rupture  velocity  Is  assumed  to  be  constant  and  subsonic.  While 
this  behavior  approximately  corresponds  to  a  particular  fracture  crit¬ 
erion  (l.e.,  the  strain  weakening  model  of  Andrews,  1976),  the  resulting 
source  model  adequately  describes  ruptures  which  grow  with  weakly 
varying  rupture  velocities.  The  rupture  velocities  determined  from 
the  modelling  will  then  represent  average  rupture  velocities. 

Our  source  models  are  kinematic  descriptions  of  the  relative  slip 
velocity  which  solve.  In  an  approximate  fashion,  the  dynamic  problem 
outlined  above.  These  "quasi-dynamic"  models  Incorporate  two  basic 
physical  considerations.  First,  because  of  the  abrupt  stress  release 
of  the  frictional  model,  the  crack  tip  Is  the  dominant  source  of  radi¬ 
ated  energy.  The  dynamic  solution  of  Kostrov  (1964)  for  rupture  growth 
with  a  constant  rupture  velocity  gives  the  self -similar,  "elliptical" 
slip  distribution,  In  which  the  slip  velocity  Is  at  a  maximum  as  the 


crack,  rip  passes  and  Chen  slows  asymptotically  to  a  constant  value. 
Second,  the  healing  of  the  rupture  starts  at  the  perimeter  of  the  rup¬ 
ture  area  as  a  result  of  the  stopping  of  the  crack  tip,  and  the  onset 
of  healing  propagates  Into  the  Interior  at  the  compresslonal  wave  velo¬ 
city.  The  healing  is  also  assumed  to  he  monotonic,  i.e.,  we  do  not 
attempt  to  Incorporate  breakout  phases  which  result  from  the  Inter¬ 
action  of  the  rupture  with  a  free  surface  CBurridge  and  Halliday,  1971) , 
and  we  presume  that  In  the  presence  of  a  finite  dynamic  frictional  stress 
the  rupture  will  heal  without  reversal  of  slip. 

These  dynamic  considerations  for  our  kinematic  models  naturally 
divide  the  slip  at  any  point  within  the  rupture  area  into  two  distinct 
parts,  which  we  will  refer  to  as  the  rupture  growth  and  the  healing. 

It  should  be  noted  that  this  healing  behavior  Is  most  directly  incorpor¬ 
ated  Into  descriptions  of  the  relative  slip  velocity,  rather  than  the 
relative  dislocation.  This  Is  particularly  opportune  considering  the 
usual  far— field  source  representation,  where  the  pulse  shapes  are  deter¬ 
mined  from  the  Integral, 

0  (x,t)  -  //  Au (£ , t-TC (x , € ) ) d£ ,  (1) 

I 

over  the  fault  surface  L  of  the  relative  slip  velocity  Au(£,t)  using 
the  travel  time,  Tc(x, C) ,  between  the  source  point  and  the  receiver, 
to  determine  the  correct  delay. 


6. 


THE  QUASI-DYNAMIC  MODEL 


g-fnomatlc  Description 

During  the  rupture  growth,  the  relative  slip  velocity  at  a  point 
£  from  the  hypocenter  is  given  by 

-  At/CtHelVv2)1*  t  >  | C |/v  (2) 

where  v  is  the  rupture  velocity  and  A  is  a  slip  velocity  which  depends 
on  the  dynamic  stress  drop,  x  ,  and  the  rupture  velocity  approximately 
as  A  =  vt^/u  (Dahlen,  1974) .  Figure  1  shows  the  resulting  slip  velocity 
as  a  function  of  time  for  a  representative  point  on  the  fault.  The 
continuation  of  the  elliptical  slip  distribution  Is  shown  as  a  dotted 
line.  Note  that  for  a  rupture  which  grows  steadily,  the  slip  velocity 
at  any  point  of  the  rupture  area  does  not  decrease  below  the  velocity 
A. 

In  order  both  to  stop  the  rupture  growth  and  to  approximate  the 
causal  healing  described  in  this  introduction,  we  multiply  this  slip 
distribution  by  the  function 

1  t<  T  (.5) 

—  a  - 

TaC5)  <  t  <  ThCC) 


hCC,t) 


(31 


7. 


where  T  CO  is  the  time  of  the  onset  of  healing  for  the  point  T  CO 

®  -  «  Q  „ 

is  the  healing  time  Ci.e.,  AuCC, T^COi  “  0),  and  h  -  T^CO  -  TSC5)  is 
the  healing  interval,  during  which  the  slip  velocity  decreases  to  zero. 

For  all  the  models  used  in  this  paper,  the  onset  of  healing, 

TgCO,  is  given  by 

V5>  -  To  -  5  -  l5-*0l/a  (4) 

where  Tq  is  the  time  from  the  nucleation  to  the  complete  healing  of  the 
rupture,  a  is  the  compresslonal  wave  velocity,  and  xq  Cthe  position  of 
the  last  point  to  heal  on  the  fault  plane)  is  a  vector  which  determines 
the  direction  and  relative  extent  of  the  asymmetry  of  the  final  rupture 
area.  The  healing  function  h(5,t)  thus  describes  a  smoothed  Cby  h) 
circular  support  function  for  the  rupture  which  is  Imploding  at  the 
compresslonal  wave  velocity.  The  interaction  of  this  function  with  the 
growing  rupture  produces  a  source  with  a  finite  rupture  area  which 
heals  into  the  interior  of  the  fault  as  desired. 

The  final  rupture  area  of  these  models  is  approximately  elliptical, 
with  eccentricity  e  <  .4,  which  may  be  varied  by  varying  xq.  In  Figure 
2,  snapshots  of  the  relative  slip  velocity  of  a  strongly  asymmetrical 
version  of  the  model  are  shown.  The  slip  velocity  has  been  smoothed 
so  that  is  may  be  readily  plotted.  The  regions  where  the  rupture  is 
healing  are  shown  darkened.  Note  how  this  rupture  is  intermediate 
between  a  circular  rupture  and  a  unilaterally  propagating  one. 
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Healing 

During  the  healing  Interval,  the  slip  velocity  at  £  is  given  by 

4u«,t)  -  AtCThC5)-t)  T.«l  <  t  <  ThC«.  C5) 

KCt2-uii/v2)'s 

For  computational  simplicity,  the  healing  interval  h  is  assumed  to  be 
constant  over  the  rupture.  This  description  of  healing  does  not  in¬ 
corporate  any  of  the  complex  diffraction  effects  (i.e.,  inhomogeneous 
wave  effects)  which  generally  occur  as  a  result  of  a  deceleration  of 
the  crack  tip  (.Madariaga,  1977;  Achenbach,  1978).  However,  the  causal 
Initiation  of  healing  which  propagates  into  the  interior  of  the  rupture 
area  at  the  compressional  wave  velocity,  is  characteristic  of  both 
in-plane  and  circularly  symmetric  numerical  fracture  models  (Madariaga, 
1976)  . 

The  assumption  that  the  P-wave  stopping  phase  Initiates  the  healing 
is  a  natural  result  of  the  frictional  model  we  are  using.  While  the 
fault  is  in  motion,  the  shear  stress  acting  across  the  rupture  surface 
is  the  dynamic  frictional  stress.  The  self -similar  solution  of  Kostrov 
(1964)  Implies  that  motion  will  continue  until  Information  (carried  by 
the  P-wave  stopping  phase)  concerning  the  finiteness  of  the  fault 
reaches  it.  The  delay  between  the  stopping  event  and  the  arrival  of 
the  stopping  phase  produces  an  overshoot  in  the  distribution  of  slip 
on  most  of  the  rupture  area;  that  is,  at  the  arrival  of  the  stopping 
phase,  the  slip  is  generally  greater  than  the  static  slip  distribution 
which  would  result  from  a  stress  drop  equal  to  the  dynamic  stress  drop. 
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Because  the  shear  stress  remains  at  the  dynamic  frictional  level  during 
healing,  the  slip  cannot  reverse  direction.  After  the  rupture  has 
healed,  the  shear  stress  across  the  fault  surface  re-adjusts  so  that 
generally  the  static  stress  drop  is  greater  than  the  dynamic  stress 
drop.  In  healing,  the  kinetic  energy  of  the  fault  motions  is  mostly 
dissipated  in  frictional  heating,  although  some  is  radiated  seism inly . 

In  stress  release  models  where  the  dynamic  frictional  stress  is  zero, 
the  slip  reverses  slightly  during  healing,  oscillating  around  the 
static  offset  whose  stress  drop  equals  the  dynamic  stress  drop.  The 
kinetic  energy  of  the  fault  motion  is  radiated  selsmicly,  damping  the 
oscillations  (Burrldge,  personal  communication) . 

In  Figure  3,  shapshots  of  the  relative  dislocation  for  the  same 
model  as  Figure  2  are  shown  with  the  healing  regions  again  darkened. 

The  rupture  area  to  the  left  of  the  healing  region  is  healed.  In  the 
topmost  graph,  the  slip  is  at  the  final  static  offset.  Note  the  slight 
cusp  to  the  right;  this  is  exactly  at  xq  and  results  from  the  uniform 
healing  behavior  for  the  whole  rupture  area:  in  Madariaga's  (1976) 
results,  the  last  point  to  heal  decelerated  more  rapidly  than  the  rest 
of  the  fault. 

Stopping  Behavior 

The  constant  healing  Interval  for  the  entire  rupture,  combined  with 
the  constant  rupture  velocity,  produces  a  perimeter  region  of  width 

A?  r  h  ctv 
v+a 


C61 
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"  across  which,  the  particle  velocity  at  the  crack,  tip  decreases  contin¬ 
uously  to  zero  as  described  by  eq.  (.51 ,  but  where  the  beginning  of  the 
Interval  Is  given  by  kl/v  rather  than  T  CCl  .  This  interaction  pro— 
duces  a  reasonable  approximation  of  a  rupture  which  stops  gradually. 
Since  h  is  a  free  parameter,  this  allows  us  to  model  ruptures  which 
stop  in  an  arbitrarily  gradual  or  sudden  fashion  by  varying  A.£,  thereby 
obtaining  body  wave  spectra  which  show  spectral  falloffs  u  with 
2  <  y  <  3  and  the  general  two  corner  frequency  envelopes  of  the  type 
discussed  by  Boatwright  C1978h)  •  The  possibility  of  modelling  this 
range  of  stopping  behaviors  represents  a  significant  advance  in  source 
modelling.  As  a  motivation  for  this  kinematic  stopping  behavior,  note 
that  for  this  region  the  function  h(?»t)  can  be  written  as 

(Th(C)-t)  -  (Th(5)-|c|/v)  (Th(g)-t)  ,  (7) 

h  h  (Th(C)-U|/v) 

where  we  may  identify  the  first  term  as  an  appropriate  decrease  of  A 
resulting  from  a  decrease  in  the  dynamic  stress  drop.  The  second  term, 
(ThC5)-t)/CT^C5)-|5|/v) ,  contains  the  time  dependence  of  the  healing  we 
have  used  for  the  interior  rupture  area. 

As  shown  in  the  series  of  snapshots  in  Figure  4,  detailing  the 
change  of  the  slip  velocity  on  this  fault  perimeter,  the  interaction  of 
the  rupture  growth  and  healing  produces  a  realistically  gradual  stopping 
at  the  fault  boundary.  As  a  result,  the  slip  distributions  of  these 
models  are  generally  more  smooth  than  those  of  constant  stress  drop 
models.  This  smoothness  is  evident  in  Figure  3,  where  the  dislocation 
at  the  healed  edge  of  the  rupture  decreases  gradually  to  zero. 
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Using  Che  technique  of  Andrews  (1974) ,  we  have  calculated  the 
static  stress  drop  for  one  of  these  models.  The  results  are  displayed 
in  Figure  S.  The  effect  of  the  smoothed  dislocation  is  to  smooth  the 
stress  drop,  weakening  the  stress  concentration  at  the  perimeter  of  the 
fault.  The  peaked  behavior  of  the  stress  drop  at  xq  is  a  result  of  our 
healing  specification,  as  mentioned  above.  A  variation  of  R  over  the 
Interior  of  the  fault  could  be  determined  such  that  the  static  stress 
drop  was  smooth  at  xq  but  this  would  have  little  effect  on  the  radiated 
pulse  shapes. 


ENERGY  FLUX  PULSE  SHAPES 


A  significant  feature  of  the  “quasi-dynamic “  source  models  is  that 
the  seismic  energy  radiated  is  consistent  with  the  theoretical  results 
of  Kostrov  (.1974)  and  Madariaga  0.976)  for  general  frictional  models  of 
rupture,  as  described  In  the  Introduction.  This  agreement  allows  us  to 
use  both  the  energy  flux  pulse  shapes  and  the  time  integrated  energy 
flux  (the  radiated  energy  per  unit  area)  of  the  body  wave  arrivals  as 
modelling  constraints.  The  energy  flux  across  a  surface  at  x,  in  a  body 
wave  travelling  with  velocity  c(x),  normally  incident 'to  the  surface, 
is  given  by 


c(x,t)  -  p(x)c(x)u2(x,t) 


(8) 


where  pGc)  is  the  density  and  uGc,t)  is  the  ground  displacement  (Bullen, 
1965,  p.  1271.  Thus,  if  the  phase  distortion  of  the  free  surface  is 
corrected  for ,  the  time  history  of  the  square  of  the  ground  velocity 
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represents  the  energy  flux  of  the  body  wave  arrival.  These  v2-plots 
were  first  used  by  Hanks  (1974)  In  an  analysis  of  the  Pacolma  Dam  re¬ 
cording,  and  more  recently,  by  Boatwright  (1978a)  In  an  analysis  of  a 
multiply  recorded  aftershock  of  the  1975  Orovllle,  California  earth¬ 
quake.  Because  of  their  non-linear  signal  enhancement,  v2-plots  pro¬ 
vide  generally  coherent,  noise-free  pulse  shapes,  suppressing  echoes  and 
low  frequency  contamination. 

The  v2-plots  can  provide  seismologists  with  particularly  useful 
waveform  constraints.  In  discussing  v 2  pulse  shapes,  we  consider  only 
the  undistorted  (elastic  whole-space}  pulse  shapes.  As  squaring  the 
ground  velocity  Is  a  non-linear  operation,  any  phase  distortion  must  be 
corrected  In  the  ground  velocity  before  squaring.  The  (undistorted) 
far -field  velocity  pulse  shapes  from  simple  sources  are  made  up  of  one 
positive  and  one  negative  pulse  with  equal  areas;  therefore  the  re¬ 
sulting  v2-plots  show  two  distinct  pulses  separated  by  an  actual  zero. 

He  will  refer  to  the  first  pulse  as  the  rupture  phase,  as  It  details 
the  growth  of  the  rupture  and  refer  to  the  second  pulse  as  the  healing 
phase. 

There  are  three  distinct  pulse  measurements  which  may  be  obtained 
from  v2-plots.  The  first  of  these  measurements,  the  width  of  the  rup¬ 
ture  phase,  provides  an  estimate  of  the  pulse  rise  time  (or  the  first 
half-cycle  of  the  velocity  trace,  t^)  and  has  been  discussed  In  detail 
by  Boatwright  (1979).  This  rise  time,  affected  by  directivity,  can 
readily  be  used  to  estimate  the  duration  of  the  rupture  growth  In  the 
direction  of  the  station  at  which  the  pulse  was  observed.  A  second 
measurement,  the  separation  of  the  healing-  and  rupture-phase  peaks, 
provides  a  pulse  duration  estimate  giving  Information  about  the  geometry 


of  Che  rupture  if  the  azimuthal  distribution  of  the  observed  pulse 
shapes  is  adequate.  Finally,  the  relative  amplitude  of  the  two  phases 
establishes  a  useful  constraint  on  the  fault  motions  *  In  general,  this 
relative  energy  content  varies  substantially  over  the  focal  sphere. 

This  is  a  direct  result  of  the  difference  in  the  particle  velocity 
behavior  behind  the  crack,  tip  and  during  healing.  The  enhanced  con¬ 
structive  Interference  of  the  rupture  phase  for  body  waves  with  phase 
velocities  on  the  fault  surface  which  approach  the  rupture  velocity 
dominates  the  v2-plots  of  shear  wave  pulse  shapes  radiated  along  the 
fault.  In  directions  near  the  normal  to  the  fault  plane  the  healing 
phase  increases  in  amplitude  and  in  these  directions  the  relative  ampli¬ 
tude  of  the  phases  is  very  sensitive  to  the  rupture  velocity  and  the 
rupture  geometry.  For  the  two  events  discussed  in  the  second  part  of 
this  paper,  this  particular  feature  of  the  waveform,  combined  with  an 
approximate  description  of  the  rupture  geometry,  is  used  to  estimate 
the  rupture  velocity  and  is  therefore  a  crucial  aspect  of  the  analysis. 

MODELLING  OF  TWO  SHUMAGIN  ISLANDS  EARTHQUAKES 

On  April  6,  1974,  two  moderate  size  (m^  “  5.8,  6.0)  earthquakes 
occurred  within  a  local  network  of  short  period  seismograph  stations 
(run  by  Lamont -Doherty  Geological  Observatory)  in  the  S hums gin  Islands, 
Alaska.  They  were  followed  by  69  recorded  aftershocks  over  the  next  two 
weeks.  Both  main  shocks  triggered  a  strong  motion  accelerograph  (SMA1) 
at  Sand  Point,  50  km  NNff  of  their  epicenters.  In  Figure  6,  we  show  a 
map  view  of  the  epicentral  area. 
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A  thorough,  analysis  of  the  sequence  la  presently  in  preparation 
(House  and  Boatwright}  which,  includes  a  full  discussion  of  the  after* 
shock,  locations,  the  sequence  b-value,  the  fault  plane  solution,  the 
far— field  waveform  analysis  and  the  source  modelling.  We  have  Included 
only  the  far— field  waveform  analysis  and  the  source  modelling  in  this 
paper,  as  an  example  of  the  application  of  the  quasi-dynamic  sources 
for  energy  flux  modelling. 


SMA1  Waveform  Anal\ 


The  SMAl  records  from  the  0153  and  0356  events  were  photographically 
enlarged,  digitized  and  instrument  corrected  using  the  technique  dis¬ 
cussed  by  Boatwright  (1978b) .  Since  the  Sand  Point  station  was  at  an 
SH  node,  the  vertical  and  horizontal  components  were  combined  to  obtain 
the  Incident  SV  pulse  shape,  using  the  free  surface  transformation. 


UsyCt) 


cos  2j  u  (t)  +  sin  j  u  Ct) . 
2cos  J  2 


(9) 


Here  j  («  33°)  Is  the  angle  of  incidence  of  the  S-wave,  u^(t)  is  the 

horizontal  component  (positive  away  from  the  source)  and  u  (T)  Is  the 

z 

vertical  component  (positive  downward) .  This  (real)  transformation  was 
derived  from  Chapter  5,  problem  5.6,  of  Akl  and  Richards  (1970).  Because 
the  transformation  Is  essentially  a  rotation  Into  the  particle  motion  of 
the  incident  SV  wave.  It  suppresses  the  evanscent  SP  arrival  expected  at 
this  range  (Anderson,  1976}  by  a  factor  of  8. 
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The  resulting  SV  acceleration  traces  for  both,  events  are  shown 
in  Figures  7  and  8,  along  with  their  respective  velocity  and  v2  traces. 
The  Integration  to  velocity  was  performed  using  a  parabolic  baseline 
technique  detailed  in  Perez  et  ad.  0-979).  Both  accelerograms  show  a 
significant  12  Hz  site  response  which  is  substantially  reduced  In  the 
integration  to  velocity. 

The  v2-plots  of  the  SV  phases  at  Sand  Point  are  remarkably  sim¬ 
ilar  in  shape  and  amplitude,  although  the  0153  pulse  shape  Is  notlcably 
more  impulsive.  This  similarity  indicates  that  the  events  probably 
share  nearly  the  same  rupture  geometry,  as  they  have  the  same  focal 
mechanism.  The  large  relative  amplitude  of  the  healing  phase  suggests 
either  that  the  rupture  propagated  towards  the  Sand  Point  station 
Cdowndlp)  or  that  the  rupture  velocity  was  slow,  about  .6  of  the  shear 
wave  velocity,  if  the  rupture 'was  approximately  circular.  These  two 
possibilities  are  considered  in  the  discussion  of  the  WWSSN  short  period 
data.  The  final  model  v2-plots  are  shown  along  with  the  data  as  dashed 
lines. 

In  Figures  9  and  10,  we  show  the  displacement  spectra  for  the  two 
events,  as  well  as  the  final  model  spectra  (dashed  lines) .  The  data 
has  been  corrected  for  attenuation  assuming  a  shear  wave  Q  of  300.  The 
site  amplification  at  10-15  Hz  shows  up  very  strongly  in  these  spectra. 
The  corner  frequencies,  marked  by  dots,  were  estimated  assuming  this 
amplification  to  be  spurious.  The  corrected  velocity  spectra  were 
Integrated  to  obtain  the  integral  of  the  squared  ground  velocity,  which 
we  will  call  Igy. 

In  order  to  estimate  the  source  dimension,  we  have  used  three 
different  measurements  of  the  SHAl  pulse  shapes  and  spectra;  i.e. ,  the 
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corner  frequency,  the  characteristic  frequency  and  the  pulse  rise  time. 
The  necessary  spectral  and  pulse  parameters  are  listed  In  Table  1. 
Following  Boatwright  0-979) ,  the  characteristic  frequency,  n^,  Is 
defined  as 


where  u^  Is  the  low  frequency  asymptote.  This  spectral  measurement 
provides  an  estimate  of  the  source  dimension  (radius  “  a)  through  the 
relation. 


K 


6 


Oil 


where  v  Is  the  rupture  velocity  and  ■  1.9  for  this  takeoff  angle. 

The  pulse  rise  time,  t^,  defined  here  to  be  the  measurable  width  of  the 
first  pulse  of  the  v2-plot  (the  rupture  phase) ,  may  also  be  used  to 
estimate  the  source  dimension,  from  the  empirical  relation  (Boatwright, 
1979) , 


t,  -  03-120  a  . 
*  16  v 


02) 


Here,  £  •  v_ sin  6  Is  the  ratio  of  the  rupture  velocity  to  be  the  phase 
c 

velocity  of  the  ray  along  the  fault  surface.  The  results  from  these  three 
source  dimension  estimates  are  shown  In  Table  2.  If  we  assume  the  ratio 
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of  the  rupture  velocity  to  the  shear  wave  velocity  to  be  bounded  as 

.55  <  2L  <  *8*  then  ve  obtain  the  first  order  estimates  of  1.0  <  a  <  1.4 

6 

km  for  the  0153  event  and  1.4  <  a  <  2.0  km  for  the  0356  event. 

WKSSN  Waveform  Analysis 

Before  fitting  for  a  particular  rupture  model,  It  Is  necessary  to 
Investigate  the  rupture  geometry  of  the  0356  event.  Ve  have  analysed 
9  short-period  P-wave  arrivals  from  6  WKSSN  stations.  The  steps  of  this 
analysis  are  shovn  In  Figure  11,  using  the  P-wave  recorded  at  GDE  as 
an  example.  The  lowest  trace  Is  the  seiaoogram  as  digitized,  and  the 
trace  above  it,  the  seismogram  after  filtering  with  a  zero-phase  band¬ 
pass  filter.  This  filter  is  made  up  from  a  triangle  smoothing  operator 
and  a  second  order  Butterworth  high  pass  filter  (corner  at  .3  Hz)  run 
forwards  and  backwards  on  the  trace.  The  third  trace  Is  the  ground 
velocity,  obtained  by  a  recursive  deconvolution  scheme  derived  from  a 
bilinear  approximation  (Kanasewlch,  1975,  p.  194)  to  the  coupled 
galvanometer-seismometer  response.  The  uppermost  trace  is  the  square  of 
the  ground  velocity.  It  1s  the  variation  of  these  v2-plots  which  we  will 
use  to  determine  the  rupture  geometry  of  the  event. 

In  Figure  12,  we  show  bow  the  pulse  shapes  vary  over  the  focal 
sphere.  The  P,  pP,  and  sP  takeoff  directions  are  plotted  relative  to 
the  fault  plane  obtained  by  House  and  Boatwright,  so  that  the  fault 
plane  Is  the  horizontal  plane  of  the  stereonet.  To  account  for  the 
pulse  shape  differences  between  the  P  and  S  body  waves,  the  phase  leaving 
as  P-waves  have  been,  corrected  to  the  appropriate  takeoff  direction 
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(open  circles)  for  S -waves  having  the  same  phase  velocity  along  the 
fault  surface,  which  requires  the  far-fleld  pulse  shapes  to  be  identi¬ 
cal  .  From  Kostrov  (1970) ,  this  corrected  takeoff  angle  Is 

9  -  sin-1  (8  sin  6  )  (.13) 

s  —  p 

a  r 

where  9  and  9  are  measured  from  the  normal  to  the  fault  plane.  Note 
p  o 

that  almost  all  the  corrected  takeoff  angles  lie  between  20°  and  30s. 
The  line  across  the  stereonet  marks  the  strike  of  the  Benloff  zone. 

In  order  to  evaluate  the  variation  of  these  v2-plots,  In  Figure  13 
we  show  the  synthetic  variation  of  the  v2-plots  over  our  model  range. 
For  these  synthetics,  we  have  used  both  circular  and  slightly  asymmet¬ 
rical  versions  of  the  quasi-dynamic  models,  fixing  the  ratio  . 

v 

In  the  left-hand  colum  of  Figure  13,  we  show  the  variation  of 
the  v2  pulse  shape  for  different  rupture  velocities  where  the  takeoff 
angle  of  the  ray  Is  at  30s  from  the  normal  to  the  fault  plane.  These 
same  synthetics  may  also  be  used  to  describe  the  variation  of  the  v2 
pulse  shapes  for  different  takeoff  angles.  For  a  circular  rupture  with 
velocity  v  ■  .66,  the  upper  and  lower  figures  approximate  the  v2  pulse 
shapes  at  35°  and  25s,  respectively,  from  the  fault  normal.  The  con¬ 
trasting  interpretations  of  these  synthetics  results  from  the  approxi¬ 
mate  similarity  of  pulse  shapes  having  similar  values  of  v_  sin  9,  for 

c 

similar  rupture  models  (Boatwright,  1979). 

In  the.  right-hand  column,  we  detail  the  pulse  shape  variations  for 
different  directions  of  asymmetrical  rupture  growth  relative  to  the 
observer.  In  these  synthetics,  the  takeoff  angle  Is  30s  from  the  fault 
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normal  and  Che  rupture  velocity  is  v  *  .66.  The  asymmetrical  models 
have  about  an  1QZ  unilateral  rupture.  If  £^($)  is  the  distance  from 
the  hypoc enter  to  the  perimeter  of  the  fault  of  the  $  direction,  the 
percent  unilateral  rupture  la  given  as 

100  •  max  5f«>)-%Ofr+iO  .  (14) 

*  5f  (*)-Kf  OH*) 

For  an  elliptical  fault  whose  hypoc enter  Is  at  one  of  the  foci,  the 
percent  unilateral  rupture  Is  simply  100  •  e,  where  e  Is  the  eccentri¬ 
city. 

In  Figure  13,  much  of  the  more  striking  variations  of  the  v2-plots 
appears  to  result  from  differences  In  the  crustal  structure  beneath  the 
stations  used.  The  Sand  Point  v2  pulse  shape  has  been  plotted  with  the 
same  time  scale.  In  order  to  show  the  relative  attenuation  present  in 
the  short  period  WWSSN  data.  In  particular,  it  is  necessary  to  point 
out  the  broadening  (perhaps  due  to  attenuation)  of  the  HKC  pulse  shapes, 
with  respect  to  the  pulse  shapes  at  nearly  the  same  takeoff  direction. 
Also  the  pulses  at  GDH  and  AAM  appear  to  have  a  crustal  reverberation 
which  is  Interfering  destructively  with  the  healing  phase  of  the  v2 
pulse  shape.  This  Interference  may  be  seen  In  the  plot  of  the  WKSSN 
analysis  In  Figure  11  as  well. 

The  slight  difference  (In  relative  amplitude  and  timing)  of  the 
depth  phase  (updip)  pulse  shapes  relative  to  the  downdip  pulse  shapes 
suggests  that  the  rupture  had  a  slight  downdip  component  of  unilateral 
rupture  (about  5Z)  and  a  rupture  velocity  of  about  v  «  .66.  Because 
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of  the  narrow  band  of  the  short  period  WWSSN  instruments  and  the 

unknown  crustal  structure  beneath  the  stations  whose  P -waves  were 

analysed,  we  can  use  these  results  only  qualitatively.  However,  It  Is 

Important  to  note  that  nearly  all  of  the  v2-plots  fall  within  the  model 

range  spanned  by  Figure  13,  from  which  we  have  determined  the  variance 

of  our  source  parameter  estimates.  As  we  have  already  estimated  the 

ratio  ja  for  these  events,  specifying  an  approximate  rupture  velocity 
v 

and  rupture  geometry  thus  determines  our  final  models. 

Final  Source  Models 

For  the  final  source  modelling  we  have  used  synthetics  generated 

by  two  circular  versions  of  the  quasi-dynamic  models.  Since  the  dis- 

_2 

placement  spectra  from  these  events  falloff  faster  than  cu  ,  we  presume 
that  the  rupture  stopped  gradually  rather  than  abruptly  (Madariaga, 
1978),  and  this  gradual  stopping  is  incorporated  into  the  models. 

Since  both  events  were  fit  with  circular  models,  the  rupture 
velocities  of  the  two  models  are  slightly  different:  v  «  .68  for  the 
0153  event  and  v  ■  .556  for  the  0356  event;  similar  results  would  have 
been  obtained  if  we  had  fixed  the  rupture  velocity  and  used  asymmetrical 
rupture  geometries.  Using  a  slightly  asymmetrical  model  and  a  rupture 
velocity  of  .66  for  the  0356  event,  lowers  the  stress  estimates  by 
=  202. 

The  source  parameters  obtained  from  our  model  fits  are  listed  In 
Table  3.  We  have  calculated  the  moment  and  radiated  seismic  energy 
using  the  formulae, 
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where  p(5  ),  p(x)  -  3.4,  2,5  gm/cm3  and  6(5  ),  8(x)  -  4.4,  2.5  km/sec 

.0  ,0 

are  the  densities  and  shear  wave  velocities  at  the  source  and  receiver, 

respectively,  i?(x,5Q)  -  52  km  is  the  geometrical  spreading  factor, 

SV 

calculated  following  Newman  (1973)  and  F  (9,$)  -  .46  is  the  radiation 
pattern  coefficient.  e^CBO8)  *.2  is  the  fractional  energy  flux,  which 
relates  the  time  Integrated  energy  flux  at  a  particular  takeoff  angle 
to  the  total  radiated  seismic  energy  (Boatwright,  1979)  . 

The  dynamic  stress  drop  also  may  be  calculated  directly  from  this 
modelling.  Since  the  quasi-dynamic  models  incorporate  the  self —similar 
slip  distribution  described  by  Kostrov  (1964) ,  Burridge  and  Willis  (1969) 
and  others,  the  slip  distribution  is  scaled  by  the  initial  relative  slip 
velocity,  A.  Thus,  any  model  fit  may  also  be  used  to  determine  the 
dynamic  stress  drop,  via  the  formula  (Boatwright,  1979), 
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synthetics,  nCx,t).  The  model  fits  give  .  uGc,t)  .  ■  .38  cm/ sec  and 

^  SCx.tl  ^ 

.34  cm/ sec  for  the  0153  and  0356  events,  respectively. 

The  results  compiled  in  Table  3  show  two  systematic  anomalies. 

For  both  events,  the  dynamic  stress  drop  is  greater  than  the  average 
static  stress  drop,  while  the  apparent  stress  is  substantially  lower 
than  x  / 4,  which  is  the  expected  value  for  frictional  ruptures  with 
v  z  .60  (Wadariaga,  1976).  However,  the  gradual  stopping  of  these 
events  may  explain  both  anomalies.  If  the  rupture  nucleated  in  a  loc¬ 
alized  region  of  high  stress,  and  grew  beyond  this  region,  the  average 
stress  drop  over  the  rupture  area  would  be  lower  than  the  initial  stress 
drop,  while  the  radiated  energy  would  be  low  due  to  the  gradual  stopping 
in  the  less  loaded  region.  Note  also  the  consistent  differences  in  the 
stress  estimates  for  the  two  events,  that  the  larger  0356  event  has  a 
stress  drop  =  70Z  of  the  0153  event. 

CONCLUSIONS 

The  first  goal  of  this  modelling  effort  is  the  determination  of 
the  rupture  velocity  of  these  two  events.  The  estimate  of  rupture  velo¬ 
city  is  critical  both  for  the  modelling,  because  of  the  trade-off  in  the 
directivity  between  rupture  velocity  and  rupture  geometry,  and  for  the 
estimates  of  stress  release  CAc  and  ,  because  they  depend  strongly  on 
the  estimate  Cor  a  priori  assumption!  of  the  rupture  velocity. 

To  illustrate  the  second  point,  note  that  if  the  rupture  velocities 
were  estimated  to  be  .98,  then  both  the  static  and  dynamic  stress  drop 
estimates  would  be  approximately  30Z  of  the  values  listed  in  Table  3. 


This  non-linear  dependence  obviously  requires  a  strong  estimate  of  the 
rupture  velocity  in  order  to  obtain  reliable  estimates  of  stress  re¬ 
lease. 

In  Figure  14,  we  have  plotted  tbe  locations  of  29  aftershocks  of 
the  tw  events,  projected  onto  the  fault  plane.  Note  the  fit  of  the 
estimated  rupture  perimeter  of  the  0153  event  to  Its  aftershock  cluster. 
The  aftershock  cluster  of  the  0356  event  Is  more  diffuse,  perhaps  be¬ 
cause  of  Its  lower  stress  drop.  While  these  aftershock  locations  cannot 
be  used  to  confirm  our  estimates  of  the  source  dimension  Cowing  to  their 
relative  uncertainty  and  the  small  dimensions  of  the  earthquakes} ,  they 
indicate  that  the  two  events  represent  spatially  distinct  ruptures  Csee 
also  Figure  6) .  This  interpretation  Is  reinforced  by  the  systematic 
differences  in  stress  released  by  the  two  events.  Because  of  the 
similarity  of  the  fault  plane  solutions  and  the  SMAl  waveforms,  the 
uncertainty  of  the  ratio  of  any  estimate  of  stress  release  is  approxi¬ 
mately  1QZ,  and  therefore  these  differences  are  significant. 

These  arguments  lead  naturally  to  the  conclusion  that  the  two 
main  shocks  represent  the  failure  of  two  distinct  asperities,  or  stress 
concentrations.  Thus,  the  extremely  high  stress  drops  are  not  directly 
Indicative  of  a  similarly  high  average  stress  over  the  region  although 
these  events  may  be  presumed  to  load  the  unruptured  part  of  the  fault 
plane.  The  stress  concentrations  are  inferred  to  be  the  result  of 
patches  of  high  strength  which  have  not  yielded  with  the  rest  of  the 
fault  plane.  Andrews  0-375}  has  called  these  stress  concentrations 
anti-dislocations,  as  they  represent  distributions  of  negative  slip 
relative  to  the  rest  of  the  fault  system.  The  gradual  stopping  charac— 
ter  of  both  events  C inferred  from  the  u  falloff  of  their  displacement 
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spectra  and  Incorporated  into  the  final  source  models)  also  supports 
this  interpretation,  as  this  stopping  behavior  Is  to  be  expected  of  a 
rupture  which  grows  beyond  the  localized  stress  concentration  where  it 
nucleated.  This  is  equivalent  to  the  "seismic  gap"  stopping  c« 

discussed  by  Hussein!  et  al .  0.976).  Using  their  formula  for  fracture 
energy. 


y  -  at  2 
o  a 

2irp 


081 


where  a  Is  the  fault  radius,  we  obtain  y  s  3  x  loiQ  ergs/cm2  for  the 
two  events.  This  extremely  high  fracture  energy  implies  that  the  fault 
zone  for  these  events  had  a  similarly  high  fracture  strength,  and  fur¬ 
ther  corroborates  the  asperity  interpretation. 
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TABLE  1  —  Spectral  Parameters 


0153  0356 


usv 

.13  ±  .04 

cm/sec 

.25  ±  .07 

cm/ sec 

*5* 

2.2  ±  .6 

cm2/sec 

3.1  ±  .8 

cm2/ sec 

1.2  Hz 

.8  Hz 

ns 

5.1  Hz 

3.5  Hz 

T% 

.24  sec 

.37  sec 
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TABLE  2  —  Rupture  Duration 


b) 


\ 


Method 

0153 

0356 

corner  frequency 

.41 

±  .1  sec 

.62 

±  .2  sec 

characteristic  frequency 

.37 

±  .1  sec 

.53 

±  .15  sec 

rise  time 

.39 

±  .08  sec 

.60 

±  .1  sec 

average 

.39 

±  .05  sec 

.59 

±  .1  sec 
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TABU  3  —  Source  Parameters 

0153  event 

radius  -  a  »  1.2  km 

moment  -  M  ■  3.5  ±  .8  *  1024  dyne-cm 
o 

static  stress  drop  -  Ac  ■  890  bars  range  600-1100  bars 
dynamic  stress  drop  -  t  •  1040.  £  350  bars 
radiated  energy  -  E  ■  8.7  ±  3.0  *  102Q  dyne-cm 
apparent  stress  -  xa  ■  160  ±  60  bars 

0356  event 

radius  -  a  *  1.65  km 

moment  -  M  ■  6.7  ±  1.5  *  1024  dyne-cm 

static  stress  drop  -  Ao  •  650  bars  range  350-800  bars 
dynamic  stress  drop  -  t&  •  78Q  ±  250  bars 
radiated  energy  -  E  ■  12.4  ±  4.0  *  1020  dyne-cm 
apparent  stress  -  t  ■ 


120  ±  50  bars 
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FIGURE  CAPTIONS 


Generalized  plot  of  the  quasi-dynamic  slip  velocity. 

T  is  the  rupture  arrival  time,  T  the  arrival  time  of 

the  first  stopping  phase  and  T  the  time  of  healing. 

h. 

The  slip  velocity  scales  with  the  velocity  A,  which  de¬ 
pends  on  the  rupture  velocity  and  dynamic  stress  drop 
approximately  as  A  s  v  Tg  . 

Snapshots  of  the  slip  velocity  distribution  of  a  quasi— 
dynamic  model,  with  time  increasing  from  the  bottom.  The 
slip  is  healing  in  the  darkened  regions.  The  slip  velocity 
has  been  smoothed  Chy  dx)  so  that  it  may  be  easily  plotted. 
Note  how  the  character  of  the  rupture  changes  from  a  cir¬ 
cular  rupture  to  a  unilateral  propagation. 

Snapshots  of  the  slip  distribution  of  the  model  shown  in 
Figure  2.  The  ship  has  healed  to  the  left  of  the  darkened 
regions.  In  the  top  figure,  the  static  slip  distribution 
Is  shown;  note  the  smoothed  slip  at  the  rupture  perimeter, 
and  the  cusp  at  xq. 
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Figure  4; 


Figure  5: 


F*Sure  6: 


Figure  7; 


Figure  8: 


Snapshots  of  the  slip  velocity  distribution  at  one  edge 
of  the  model,  detailing  the  stopping  behavior.  The 
healing  regions  have  not  been  darkened.  The  stopping 
phase  can  be  seen  as  a  discontinuity  in  slope  at  times 
t-  .16  and  t  -  .18.  By  t  ■  .22  the  perimenter  has  healed 
completely. 

Distribution  of  static  stress  drop  and  final  slip  for  one 
of  the  models.  Note  how  the  smoothed  distribution  of 
slip  naturally  weakens  the  stress  concentration  at  the 
perimeter  of  the  rupture.  The  sharp  stress  node  is  an 
unphysical  artifact  of  the  healing  description. 

Map  of  the  Aleutian  arc  near  the  Shumagin  Islands,  showing 
the  epicentral  area  of  the  events  to  be  discussed  and  the 
stations  used  to  locate  the  events  and  their  aftershocks. 
The  inset  shows  the  epicenters  of  the  two  events  along 
with  the  29  located  aftershocks.  Note  the  clear  grouping 
of  the  aftershocks  into  two  distinct  clusters. 

SV  acceleration,  velocity  and  v2  pulse  shapes  for  the  0153 
event.  The  dashed  line  in  the  v2-plot  is  the  synthetic 
pulse  shape  from  the  final  model. 

SV  acceleration,  velocity  and  v2  pulse  shape  for  the  0356 
event.  The  dashed  line  in  the  v2-plot  is  the  synthetic 
pulse  shape  from  the  final  model. 
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Figure  9: 


Figure  10: 


Figure  11; 


Figure  12; 


Figure  L3 : 


Displacement  amplitude  spectrum  for  the  0153  event.  The 
dot  marks  the  measured  corner  frequency.  The  spectrum 
has  been  corrected  for  attenuation.  The  spectral  ampli- 
.  fication  at  10-15  Hz  is  a  site  response.  The  dashed 
line  is  the  synthetic  spectrum. 

Displacement  amplitude  spectrum  for  the  0356  event.  See 
expanatlon  for  Figure  9. 

Detail  of  short-period  WWSSN  analysis.  The  lowermost 
trace  is  the  seismogram  as  digitized  with  the  bandpassed 
seismogram  above  it.  The  next  traces  are  the  deconvolved 
velocity  and  finally  the  v2-plot. 

Variation  of  the  short-period  WWSSN  v2  pulse  shapes  from 
the  0356  event  over  the  focal  sphere.  The  takeoff  angles 
of  the  phases  have  been  rotated  so  that  the  fault  plane 
is  the  plane  of  the  stereonet.  The  phases  which  took 
off  as  P-waves  are  corrected  (solid  lines)  to  the  equi¬ 
valent  takeoff  angles  for  S-waves.  The  two  arcs  are  at 
20°  and  30°  from  the  fault  normal.  The  Sand  Point  SMA1 
v2-plot  is  plotted  at  the  same  time  scale  for  reference. 

The  variation  of  the  v2  pulse  shapes  over  the  model  range. 
The  left-hand  column  shows  the  variation  of  the  pulse 
shapes  with  the  variation  of  rupture  velocity  or  takeoff 
angle  (see  text).,  while  the  right-band  column  shows  the 
variation  resulting  from  a  slightly  asymmetrical  rupture 


geometry. 
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Figure  14:  Plot  of  the  aftershock  locations  projected  onto  the 

fault  plane.  The  dashed  circles  are  the  rupture  areas 
determined  from  the  waveform  modelling.  The  uppermost 
(smaller)  event  is  the  0153  main  shock,  while  the  larger 
event  is  the  0356  main  shock.  The  largest  intervening 
aftershock  (m^  •  4.7)  occurred  at  0227  and  is  just  to  the 
left  of  the  0153  rupture  zone. 
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INVESTIGATION  OF  TWO  HIGH  STRESS-DROP  EARTHQUAKES 
IN  THE  SHUMAGIN  SEISMIC  GAP,  ALASKA1 

Leigh  House  and  John  Boatwright 

Lcunont.-Voh.eAty  Geo  tog  teat  ObtenvatoKy  and 
department  o  j  Geotogtcat  Science*  otf  Cotumbta  Untverutty 
Patt&ade*,  New  Yonk  10964 

ABSTRACT 

Two  moderate  size  (m^  ■  5.8,  6.0)  earthquakes  occurred  within  a  local  net¬ 
work  of  short-period  seismograph  stations  in  the  Shumagin  Islands,  Alaska,  on 
April  6,  1974.  They  were  followed  by  69  aftershocks  recorded  over  the  next  two 
weeks.  Both  mainshocks  triggered  a  strong-motion  accelerograph  (SMA)  at  Sand 
Point,  50  km  NNW  of  their  epicenters. 

High  quality  locations  obtained  from  local  network  arrivals  for  the  main- 
shocks  and  29  aftershocks  plot  at  depths  between  35  km  and  45  km  and  define  a 
plane  dipping  about  30°  to  the  NW.  A  nearly  pure-thrust  focal  mechanism  for 
the  larger  (m^  *  6.0)  earthquake  was  obtained  from  long-period  data.  The  fault 
plane  dips  30s  in  the  direction  N  16°W.  This  sequence  was  located  along  the 
dipping  seismic  zone  beneath  the  eastern  Aleutians  and  was  presumably  related 
to  underthrusting  of  the  Pacific  plate  beneath  North  America. 

We  obtained  estimates  of  the  source  parameters  of  these  earthquakes  from 
analysis  of  SMA  data  and  WWSSN  short  period  data.  WWSSN  data  Indicates  that 
the  earthquakes  ruptured  a  nearly  circular  zone.  Modelling  of  the  SMA  records 
with  a  quasi-dynamic  model  [Boatwright,  1979]  provides  the  following  source 
parameter  estimates  for  the  m^  ■  5.8  and  6.0  earthquakes  respectively:  moments. 
Mo,  3.6  and  6.6,  x  102**  dyne-cm  and  stress  drops:  650  and  540  bars.  A  high 
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frequency  spectral  fall-off  of  <o~3  suggests  that  the  ruptures  stopped  gradually. 

The  Shumagln  Islands  region  Is  believed  to  have  a  high  potential  for  a 
future  large  earthquake  [Kelleher,  1970;  Sykes,  1971;  Kelleher  e£  al. ,  1973]. 

The  location  of  this  earthquake  sequence  at  the  deepest  part  of  the  rupture 
zone  of  the  1938  earthquake  (Ms  ■  8.7)  (major  earthquake  ruptures  of  ten  Ini¬ 
tiate  at  depth  and  propagate  updip)  and  the  high  stress-drops  of  the  shocks  in 
1974  may  indicate  considerable  accumulation  of  stress  prior  to  a  major  earth¬ 
quake  in  the  Shumagln  Islands  region. 

INTRODUCTION 

Two  moderate  size  (m^  “  5.8,  6.0)  earthquakes  and  their  aftershocks,  which 
occurred  within  the  Shumagln  Islands  Seismic  Network,  Alaska,  have  produced  a 
unique  data  set  for  a  detailed  study  of  the  tectonics  at  depth  in  an  area  which 
has  been  identified  as  a  seismic  gap.  In  order  to  fully  describe  this  earth¬ 
quake  sequence,  this  study  Integrates  the  locations  and  magnitudes  from  the 
seismic  network,  the  strong  motion  accelerograph  recordings  of  both  events,  and 
short-period  WWSSN  data  from  the  larger  event. 

The  April,  1974  sequence  began  at  0153  hours  on  April  6,  1974,  with  a  mod¬ 
erate  size  (m^  "  5.8)  mainshock.  This  was  immediately  followed  by  aftershocks, 
including  one  of  magnitude  4.3  at  0227  hours.  A  second,  and  larger  (m^  »  6.0), 
mainshock  occurred  at  0356  hours.  Over -the  next  two  weeks  there  were  nearly  70 
aftershocks,  three  of  which  had  body  wave  magnitudes  greater  than  4.  Considering 
only  the  teleseismic  locations,  the  sequence  resembles  a  swarm  as  the  mainshocks 
are  so  close  in  both  location  and  magnitude.  Using  the  local  network  locations, 
however,  it  is  clear  that  there  were  two  distinct  mainshocks  which  had  separate 
rupture  areas.  The  range  of  hypocentral  depths  for  these  events  is  between  37 


and  43  km. 
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The  Lamont-Doherty  Geological  Observatory  has  operated  a  network  of  vertical 
short-period,  radio-telemetered  stations  in  the  Shumagin  Islands  region  of  Alaska 
since  July  1973.  As  originally  installed,  the  network  consisted  of  high-remote 
stations  which  telemetered  their  data  to  a  central  recording  site  at  Sand  Point 
(SDP,  see  Figure  1).  Severe  environmental  factors  resulted  in  numerous  station 
outages  over  the  winter  season.  By  April,  1974,  when  the  sequence  occurred,  only 
four  of  the  remote  stations  were  operating. 

Both  mainshocks  triggered  a  strong  motion  accelerograph  (SMA-1)  located  at 
Sand  Point.  These  recordings,  complemented  by  a  qualitative  analysis  of  the 
rupture  geometry  using  WWSSN  short-period  data,  are  modelled  using  the  techniques 
and  models  discussed  by  Boatwright  [1979],  which  focus  on  the  determination  of 
the  rupture  velocity,  and  which  provide  more  reliable  estimates  of  stress  re¬ 
lease  . 

Archambeau  [1978]  made  a  survey  of  earthquakes  along  the  Aleutian  arc  using 
m^  to  Ms  ratios  and  obtained  high  stress  drops  for  these  two  mainshocks  (Act  = 

500  bars).  Our  study  confirms  his  identification  of  these  events  as  high-stress 
drop  earthquakes,  and  in  addition,  provides  considerable  Insight  into  the  charac¬ 
ter  of  the  deformation  at  this  depth  in  the  Benloff  zone. 

Analysis  of  Shumagin  Network  Data 

During  the  two  weeks  following  the  mainshocks  of  April  6,  1974,  more  than 
70  aftershocks  occurred  which  were  large  enough  to  be  recorded  by  the  station  at 
Squaw  Harbor  (SQH,  see  Figure  1).  Of  these,  4  were  large  enough  (m^  >  4.0)  to 
be  recorded  teleseismically .  Seismic  stations  which  were  operating  at  the  time 
of  these  earthquakes  were:  SDP,  PW,  SQH,  SGB,  and  CNB  (see  Figure  1).  CNB  was 
working  only  intermit tantly  and  recorded  only  about  1/2  of  the  sequence.  Unfor¬ 
tunately,  the  station  nearest  the  sequence,  NG1,  did  not  record  any  of  the  sequence. 


Our  earthquake  location  procedure  consisted  of  reading  arrival  times  from 
magnetic  tape  playbacks  (with  a  precision  of  0.2s)  and  using  a  variation  of  the 
HYPO  71  program  [Lee  and  Lahr,  1972]  to  obtain  hypocentral  coordinates  from  the 


arrival  times.  We  had  no  difficulty  making  P  arrival  time  picks,  as  the  P 
phase  was  generally  quite  impulsive.  Since  only  vertical  seismometers  are  in¬ 
stalled  at  the  remote  stations,  making  reliable  picks  for  S  arrivals  was  quite 
difficult.  In  general,  we  used  only  S  arrival  times  from  the  SDP  station, 
which  has  2  horizontal  short  period  seismometers,  in  addition  to  a  vertical. 

In  order  to  reduce  the  magnitude  of  the  overall  station  residuals,  we 
applied  station  corrections  to  arrival  times.  We  averaged  the  station  resi¬ 
duals  from  the  two  main  shocks  and  the  first  large  aftershock  (April  6,  0227 
hrs,  n»k  -  4.3),  and  applied  the  corrections  to  arrivals  from  the  whole  suite 
of  events.  The  largest  residual,  -0.1s,  was  at  station  SQH;  the  rest  were 
+  0.05s  or  less. 

By  obtaining  earthquake  locations  both  with  and  without  arrival  times  from 
CNB,  we  determined  that,  although  there  is  no  significant  bias  of  locations 
which  don’t  use  CNB  arrival  times  as  compared  with  those  which  do,  there  is  an 
Increase  in  scatter  of  the  locations  which  don't  use  CNB  arrival  times  recorded 
at  CNB.  This  restriction  reduced  the  number  of  aftershocks  we  located  to  27. 

We  feel  that  these  locations  are  the  most  reliable  of  those  from  the  whole  se¬ 
quence  and  estimate  hypocentral  location  errors  to  be  less  than  5  km  in  an 
absolute  sense,  and  in  a  relative  sense,  less  than  2  km  for  our  "A"  quality 
solutions. 

Locations  of  the  April  6,  1974  mainshocks  and  well  located  aftershocks  are 

i 

plotted  in  map  view  in  Figure  1.  Two  open  circles  within  the  box  are  network 
locations  of  the  two  mainshocks.  The  PDE  (teleseismic)  location  for  these  events 
is  plotted  as  the  open  circle  about  20  km  to  the  NNW  of  their  network  locations. 


r 
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Network  locations  for  the  whole  sequence  are  plotted  In  the  inset.  Main  shock 
locations  are  circled;  the  first,  at  0153  hours  on  April  6,  is  trenchward  of 
the  second,  which  occurred  two  hours  later,  at  0356.  Note  the  separation  of 
aftershocks  into  two  adjacent,  but  distinct  groups. 

The  location  of  the  cross  section  in  Figure  2  is  indicated  in  Figure  1 
by  the  line  which  terminates  just  SE  of  the  Aleutian  Trench.  The  cross  section 
extends  slightly  beyond  the  northern  edge  of  Figure  1. 

Four  years  of  Shumagin  Seismic  Network  data  are  plotted  in  cross  section 
views  in  Figure  2.  This  figure  illustrates  the  well-defined  10  km  thin  Benioff 
zone  which  exists  within  the  centred  portion  of  the  network.  The  April  6  se¬ 
quence  occurred  within  the  area  of  the  box,  at  a  depth  of  40  km.  Davies  and 
House  [1979J  noted  that  if  the  seismicity  below  about  40  km  occurs  near  the 
upper  surface  of  the  descending  slab,  then  there  must  be  a  bend  in  the  slab  at 
about  40  km.  This  is  necessary  because  the  seismic  zone  below  40  km  dips  at 
about  30® ,  whereas,  between  the  trench  and  40  km  depth,  the  slab  dip  is  about 
15°.  Thus,  the  April  6  sequence  occurred  very  near  this  bend  in  the  slab. 

The  inset  in  Figure  2  is  an  enlarged  cross  section  of  the  area  of  the  April  6, 
1974  sequence.  The  first  main  shock,  at  0153  hours,  is  shallower,  and  trench- 
ward  of,  the  2nd  main  shock.  Dashed  lines  represent  the  rupture  dimensions  as 
estimated  from  the  source  modelling.  The  rupture  zones  are  plotted  with  a  dip 
of  30®,  which  is  the  clip  of  the  fault  planes  we  infer  for  these  events.  Sym¬ 
bol  size  is  scaled  to  magnitude,  the  largest  event  is  the  second  mainshock  (m^  “ 
6.0),  the  smallest  event  plotted  has  a  magnitude  of  about  1.5.  Symbol  filling 
indicates  quality  of  the  location.  Filled  symbols  are  used  for  the  best  quality 
locations  ("A"  quality),  1/2  filled  for  the  next  best  ("B"  quality)  and  open 
symbols  for  worst  quality  ("C").  Relative  hypocentral  errors  are  about  1.5  km 
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for  the  "A"  quality  locations,  2  km  for  "B"  quality  and  3.5  km  for  "C"  quality. 
Errors  in  "A"  quality  locations  are  about  the  same  magnitude  as  the  rupture 
radius  of  the  first  malnshock  (1.3  km),  those  for  "B"  quality  about  the  same 
magnitude  as  the  rupture  radius  of  the  second  mainshock  (1.8  km). 

Figure  3  is  an  Inclined  cross  section  view  of  the  April,  1974  sequence. 

The  projection  plane  is  oriented  parallel  to  the  fault  planes  of  the  mainshocks, 
and  the  "view"  of  this  figure  is  upwards  at  the  fault  plane.  In  this  figure 
the  first  main  shock  (0153  hours)  is  above  the  second  (0356  hours).  The  dashed 
circles  represent  the  rupture  areas  of  the  two  mainshocks  as  inferred  from  the 
SMA  data.  Symbol  size  and  filling  represent  magnitude  and  quality  of  the 
solutions  as  in  the  previous  figure.  Note  that  the  aftershocks  cluster  about 
the  two  main  shocks,  and,  in  general,  are  Indicative  of  two  distinct  rupture 
areas.  There  is,  however,  some  scatter,  which  is  probably  partly  the  result 
of  location  errors,  but  which  may  also  suggest  that  deformation  during  this 
sequence  extended  beyond  the  immediate  rupture  zones  of  the  mainshocks.  The 
first  large  aftershock  (m^  *4.3)  occurred  at  0227  hours,  on  April  6  (just 
1/2  hour  after  the  first  malnshock).  It  is  plotted  as  the  large,  filled  sym¬ 
bol  to  the  left  of,  and  slightly  below  the  first  mainshock,  as  shown  in  Figure  3. 
The  second  large  aftershock  (m^  ■  4.1)  occurred  at  0509  hours  on  April  6,  and 

plots  as  the  solid  square  to  the  left  of,  and  slightly  below,  the  2nd  mainshock. 

\ 

Thus,  the  locations  of  these  larger  aftershocks  also  suggests  that  deformation 
during  the  sequence  extended  beyond  the  immediate  malnshock  rupture  areas.  The 
scatter  in  aftershocks  plotted  in  the  inset  of  Figure  2  is  consistent  with  this 
concept . 

Focal  Mechanism  of  the  Main  Shocks 

We  obtained  the  focal  mechanism  of  the  second  mainshock  (m.  -  6.0),  shown 
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in  Figure  4,  primarily  from  long  period  arrivals  at  WWSSN  stations.  Both  S  wave 
polarizations  and  P  wave  first  motions  were  used.  There  is  only  one  inconsis¬ 
tency  in  the  first  motions,  and  that  is  a  less  reliable  pick.  This  mechanism 
is  consistent  with  the  local  network  short  period  first  motions,  which  are  the 
square  symbols  in  Figure  4. 

Since  error  in  determination  of  Che  mechanism  would  produce  error  in  the 
source  parameter  estimates  (see  below) ,  we  wanted  to  extract  the  maximum  constraint 
possible  from  the  data  and  quantify  error  in  the  focal  mechanism.  We  selected 
S  wave  polarizations  from  9  stations  which  had  clear  S  arrivals,  and  used  the 
S  arrival  from  the  SDP  SMA,  as  well,  and  obtained  a  mechanism  which  produced  a 
minimum  S  wave  polarization  residual.  We  also  used  this  information  to  obtain  1- 
standard  deviation  error  estimates.  Focal  parameters  and  error  limits  are:  strike, 
254 #  +  15°,  dip  30°  +  5°,  rake  (of  slip  vector)  90°  +  15°.  These  errors  are  indi¬ 
cated  on  Figure  4,  as  are  the  associated  errors  in  the  T  and  P  axes.  First  motion 
and  S  wave  polarization  data  are  fewer,  but  Identical,  for  the  first  malnshock 
(m^  =  5.8). 

Since  the  network  located  seismic  zone,  and  the  general  distribution  of  after¬ 
shocks  (see  Figure  2)  are  very  nearly  parallel  to  the  NE-SW  striking,  30°  dipping 
nodal  plane,  we  prefer  this  to  be  the  fault  plane  of  these  events.  The  SMA  wave 
forms  at  SDP  are  also  consistent  with  this  choice  of  fault  plane.  Thus,  these 
earthquakes  were  shallow-angle  underthrusts  of  Pacific  lithosphere  beneath  North 
America.  Shallow  angle  under  thrusting,  with  mechanisms  similar  to  this  one  has 
been  observed  along  the  Aleutians,  by  investigators  such  as  Stauder  [1968]  and 
Bollinger  and  Stauder  [1966] .  Since  their  fault  planes  are  parallel  to  the 
Benloff  zone,  as  shown  by  the  network  locations,  below  40  km,  it  appears  that 
the  Pacific  plate  has  already  made  the  bend  to  a  steeper  dipping  geometry  by  the 
time  It  has  reached  the  40  km  depth  of  these  earthquakes. 
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SMA.-1  Waveform  Analysis 

The  SMA-1  records  from  the  0153  and  0356  events  were  photographically  en¬ 
larged,  digitized  and  instrument  corrected  using  the  techniques  discussed  in 
Boatwright  (1978).  Since  the  Sand  Point  station  was  at  an  SH  mode,  the  ver¬ 
tical  and  horizontal  components  were  combined  to  obtain  the  incident  SV  pulse 
shape,  using  the  free  surface  transformation, 

u  (t)  -  u  (t)  +  sin  j  u  (t).  (1) 

sv  2  cos j  x  J  z 

Here  j  is  the  angle  of  incidence  of  the  s-wave,  u^(t)  is  the  horizontal  compo¬ 
nent  (positive  away  from  the  source)  and  u  (t)  is  the  vertical  component  (posi- 

z 

tlve  downward).  This  (real)  transformation  was  derived  from  Chapter  5,  problem 

5.6,  of  Akl  and  Richards  (1979).  The  resulting  SV  acceleration  traces  for  both 

events  are  shown  in  Figures  6  and  7,  along  with  their  respective  velocity  and 
2 

v  traces.  The  integration  to  velocity  was  performed  using  a  parabolic  base¬ 
line  technique  detailed  in  Perez,  Husid  and  Espinosa  (1979).  Both  accelerograms 
show  a  significant  12  hz  site  response  which  is  substantially  reduced  in  the 

integration  to  velocity. 

2 

The  v  -plots  detail  the  energy  flux  of  the  body  wave  arrival.  Their  non¬ 
linear  signal  enhancement  make  them  a  strong  tool  for  seismic  source  studies. 

2 

The  v  -plots  of  the  sv  phases  at  Sand  Point  are  remarkably  similar  in  shape 
and  amplitude,  although  the  0153  pulse  shape  is  noticably  more  impulsive.  This 
similarity  indicates  that  the  events  probably  share  nearly  the  same  rupture 
geometry  as  they  have  same  focal  mechanism.  The  large  relative  amplitude  of 
the  healing  phase  suggests  either  that  the  rupture  propagated  towards  the  Sand 
Point  station  (downdip)  or  that  the  rupture  velocity  was  slow,  about  .6  of  the 
shear  wave  velocity,  if  the  rupture  was  approximately  circular  (see  Boatwright, 
1979).  These  two  possibilities  are  considered  in  the  discussion  of  the  WWNSS 
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short  period  data.  The  final  model  v  -plots  are  shown  along  with  the  data  as 
dashed  lines. 

In  Figures  8  and  9,  we  show  the  displacement  spectra  for  the  two  events, 
as  well  as  the  final  model  spectra  (dashed  lines).  The  data  has  been  correc¬ 
ted  for  attenuation  assuming  a  shear  wave  Q  of  300.  The  site  amplification 
at  10-15  hz  shows  up  very  strongly  in  these  spectra.  The  corner  frequencies, 
marked  by  dots,  were  estimated  assuming  this  amplification  to  be  spurious. 

The  corrected  velocity  spectra  were  integrated  to  obtain  the  integral  of  the 
squared  ground  velocity  (say  ^ ) . 

The  spectral  and  pulse  shape  parameters  for  the  two  events  necessary  for 
the  estimation  of  the  rupture  dimensions  are  compiled  in  Table  1.  Following 
Boatwright;  (1978a),  the  characteristic  frequency,  n g  »  is  defined  as 


where  u  is  the  low  frequency  asymptote.  This  spectral  measurement  provides 
f* 

an  estimate  of  the  source  dimension  (radius  -  a)  through  the  relation. 


-(f) 


(3) 


where  v  is  the  rupture  velocity  and  tc 


S 


69  for  this  takeoff  angle.  The  rise 


time,  t^,  to  be  the  measurable  width  of  the  first  pulse  of  the  v  -plot 

(the  rupture  phase),  may  also  be  used  to  estimate  the  source  dimension,  from  the 

empirical  relation  (Boatwright,  1979), 


-  ^(f)-  <*> 

y 

Here  4  *  ~  sin  6  is  the  ratio  of  the  rupture  velocity  to  the  phase  velocity  of 
the  ray  along  the  fault  surface.  The  results  from  these  three  source  dimension 


estimates  are  shown  in  Table  2 


If  we  assume  the  ratio  of  the  rupture  velocity. 
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to  the  shear  velocity  to  be  bounded  as  .  55  <  ^-  <  .8,  then  we  obtain  the  first 
order  estimates  of  1.0  <  a  <  1.4  km  for  the  0153  event  and  1.4  <  a  <  2.0  for 
the  0356  event. 

WWNSS  Waveform  Analysis 

Before  fitting  for  a  particular  rupture  model,  it  is  necessary  to  inves¬ 
tigate  the  rupture  geometry  of  the  0356  event.  We  have  analyzed  9  short  period 
P-wave  arrivals  from  6  WWNSS  stations.  The  steps  of  this  analysis  are  shown 
in  Figure  w,  using  the  P-wave  recorded  at  GDH  as  an  example.  The  lowest  trace 
is  the  seismogram  as  digitized,  and  the  trace  above  it,  the  seismogram  after 
filtering  with  a  zero-phase  bandpass  filter.  This  filter  is  made  up  from  a  tri¬ 
angle  smoothing  operator  and  a  second  order  Butterworth  high  pass  filter  (corner 
at  .3  hz)  run  forwards  and  backwards  on  the  trace.  The  third  trace  is  the 
ground  velocity,  obtained  by  a  recursive  deconvolution  scheme  derived  from  a 

bilinear  approximation  to  the  coupled  galvanometer-seismometer  response.  The 

2 

uppermost  trace  is  the  square  of  the  ground  velocity.  It  is  these  v  -plots 

which  we  will  use  to  constrain  the  rupture  geometry  of  the  event. 

2 

In  Figure  11  we  show  the  variation  of  the  v  -plots  over  the  focal  sphere. 

The  P,  pP  and  sP  takeoff  angles  are  shown  relative  to  the  fault  plane,  rotated 

2 

into  the  plane  of  the  figure,  together  with  their  respective  v  pulse  shapes. 

To  account  for  the  pulse  shape  differences  between  the  P  and  S  body  waves,  the 
phases  leaving  as  P-waves  have  been  corrected  to  the  appropriate  take-off  angle 
(open  circles)  for  an  S-wave  having  the  same  phase  velocity  along  the  fault 
surface,  which  requires  the  far-field  pulse  shapes  to  be  identical.  From 
Kostrov  (1970),  this  corrected  take-off  angle  is 


where  8.  and  0  are  measured  from  the  normal  to  the  fault  plane.  Note  that 
S  a 
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almost  all  the  corrected  takeoff  angles  lie  between  20®  and  30®.  The  line  across 

the  stereonet  marks  the  strike  of  the  Benioff  zone. 

2 

In  order  to  evaluate  the  variation  of  these  v  -plots,  in  Figure  12  we  show 

2 

the  synthetic  variations  of  the  v  -plots  over  our  model  range.  For  these  syn¬ 
thetics,  we  have  used  both  circular  and  slightly  asymmetrical  versions  of  the 
quasi-dynamic  models  discussed  in  Boatwright  (1979),  fixing  the  ratio  —  . 

2 

In  the  left-hand  column  of  Figure  12,  we  show  the  variation  of  the  v  pulse 

shape  for  different  rupture  velocities  where  the  takeoff  angle  of  the  ray  is 

at  30°  from  the  normal  to  the  fault  plane.  These  same  synthetics  may  also  be 

2 

used  to  describe  the  variation  of  the  v  pulse  shapes  for  different  takeoff 

angles.  For  a  circular  rupture  with  rupture  velocity  v  ».6_  ,  the  upper  and 

2 

lower  figures  approximate  the  v  pulse  shapes  at  37°  and  24°,  respectively, 
from  the  fault  normal.  The  contrasting  interpretations  of  these  synthetics  re¬ 
sult  from  the  approximate  similarity  of  pulse  shapes  having  similar  values  of 

y 

—  sin  9,  for  similar  rupture  models. 

In  the  right  hand  column,  we  detail  the  pulse  shape  variations  for  differ¬ 
ent  directions  of  asymmetrical  rupture  growth  relative  to  the  observer.  In 


these  synthetics,  the  takeoff  angle  is  30°  from  the  fault  normal  and  the  rupture 
velocity  is  v  *  . 6  .  The  asymmetrical  models  have  about  an  11Z  unilateral  rup- 

9 

ture:  if  (<fr)  is  the  distance  from  the  hypocenter  to  the  perimeter  of  the 
fault  in  the  $  direction,  the  percent  unilateral  rupture  is  given  as 


100  max 


€f(<fr)  T  Cf(4>  +  *) 
+  ") 


Note  how  the  asymmetrical  models  produce  a  strong  variation  of  the  peak  sepa- 
2 

ration  in  the  v  pulse  shapes  over  the  focal  sphere,  as  well  as  varying  the  rela¬ 
tive  amplitude  of  the  peaks.  This  difference  will  also  help  resolve  the  trade- 


L2 


off  between  rupture  velocity  and  rupture  geometry. 

2 

In  Figure  11,  much  of  the  more  striking  variations  of  the  v  -plots  appears 

to  result  from  differences  in  the  crustal  structure  beneath  the  stations  used. 

2 

The  Sand  Point  v  pulse  shape  has  been  plotted  with  the  same  time  scale,  in 

order  to  show  the  relative  attenuation  present  in  the  short  period  WWNSS  data. 

In  particular,  it  is  necessary  to  point  out  the  broadening  (perhaps  due  to 

attenuation)  of  the  HKC  pulse  shapes,  with  respect  to  the  pulse  shapes  at  nearly 

the  same  takeoff  angles.  Also  the  pulses  at  GOH  and  AAM  appear  to  have  a  crustal 

reverberation  which  is  interfering  destructively  with  the  healing  phase  of  the 
2 

v  pulse  shape.  This  interference  may  be  seen  in  the  plots  of  the  WWNSS  analy¬ 
sis  in  Figure  10  as  well. 

The  slight  differences  (in  relative  amplitude  and  timing)  of  the  depth 

phase  (updlp)  pulse  shapes  relative  to  the  downdip  pulse  shapes  suggests  that 

the  rupture  had  a  slight  downdip  component  of  unilateral  rupture  velocity  (about 

5Z)  and  a  rupture  velocity  of  about  v  -  . 6  .  Because  of  the  narrow  band  of  the 

short  period  WWNSS  instruments  and  the  (unknown)  crustal  structure  beneath  the 

stations  whose  P-waves  were  analyzed,  we  can  only  use  these  results  qualitatively 

2 

However,  it  is  important  to  note  that  nearly  all  of  the  v  -plots  fall  within  the 
model  range  spanned  by  Figure  12,  from  which  we  have  determined  the  variance  of 
our  source  parameter  estimates.  As  we  have  estimated  the  ratio  ^  for  these 
events,  specifying  an  approximate  rupture  velocity  and  rupture  geometry  then 
determines  our  final  models. 


Final  Source  Models 

For  the  final  source  modelling  we  have  used  two  circular  versions  of  the 
quasi-dynamic  models.  These  models  have  the  dynamically  feasible  "elliptical" 


or  self-similar  slip  distribution  during  rupture  growth,  and  causal  healing. 


Figure  13  shows  half-view  snapshots  of  the  slip  velocity  for  this  circular  model. 
Since  the  displacement  spectra  from  these  events  falloff  faster  than  ui-2  ,  we 
presume  that  the  ruptures  stopped  gradually  rather  than  abruptly,  [Madariaga, 
1978],  and  this  gradual  stopping  is  incorporated  into  the  models. 

Since  both  events  were  fit  with  circular  models,  the  rupture  velocities  of 
the  two  models  are  slightly  different:  v  *  .6  for  the  22  0153  event  and 
v  *  .55 Q  for  the  0356  event;  similar  results  would  have  been  obtained  if  we 
had  fixed  the  rupture  velocity  using  a  slightly  asymmetrical  model  for  the  0356 
event,  and  a  rupture  velocity  of  .6  lowers  the  stress  estimates  by  20%. 

The  source  parameters  obtained  from  our  model  fits  are  listed  in  Table  3.  We 


have  calculated  the  moment  and  radiated  seismic  energy  using  the  formulae, 

,  5  ,  ,  i?(x,£  )7 

m  -  Air  p(£  r  6  (*)T  *(*)*  8  (xV  ~ — 2 —  u  (j 

°  ~°  '  ^(e.^e  6 


p  (x)  3  (x)  8 
%  (*>5 


where  p(5  ),  p  (x)  ■  3.4,  2.5  gm/cm3  and  3  (£  ),  3  (x)  ■  A. 4,  2.5  km/sec  are  the 
o  ~o 

densities  and  shear  wave  velocities  at  the  source  and  receiver,  respectively, 
A(x,€0)  -=  52  km  is  the  geometrical  spreading  factor  calculated  following  New¬ 
man  [1973]  and  Fey(0,^)  *  .46  is  the  radiation  pattern  factor.  $  (30°)  .2  is 

the  fractional  energy  flux,  which  relates  the  time  integrated  energy  flux  at  a 
particular  takeoff  angle  to  the  total  radiated  seismic  energy  [Boatwright,  1979]. 

The  effective  stress  or  the  dynamic  stress  drop,  [Brune,  1970],  also  may  be 
calculated  directly  from  this  modelling.  Since  the  quasi-dynamic  models  incor¬ 
porate  the  self-similar  slip  distribution  described  by  Hoshov  [1964],  Burridge 
and  Willis  [1969]  and  others,  the  slip  23  distribution  is  scaled  by  the  Initial 
relative  slip  velocity,  A(o).  This  slip  velocity  is  related  to  the  dynamic  stress 


by  the  approximate  relation; 


(9) 


good  for  subsonic  rupture  velocities.  Thus  any  model  fit  also  determines  effec¬ 
tive  stress  by  the  formula; 


i.  1  i  i 

p(eo)2  B(eo)2  P (x) 2  3  (x)2 


,x) 

-0  - 


u(x,t) 

n(x,  t)6 


(10) 


where 


0(x,t). 


u(x,t) 

T -  is  obtained  by  scaling  the  data,  u  (x,t)  to  the  synthtics, 

fl(x,t) 

u(x,t) 


The  model  fits  give  - -  *  .38  and  .34  for  the  0153  and  0356  events 


0(x,t) 


respectively. 

The  results  compiled  in  Table  3  show  two  systematic  anomalies.  For  both 
events,  the  effective  stress  is  greater  than  the  average  stress  drop,  while  the 
apparent  stress  is  substantially  lower  than  Tg/^»  which  is  the  expected  value 
for  frictional  ruptures  with  v  :  .6  .  However,  the  gradual  sopping  of  these 

events  may  explain  both  anomalies.  If  the  rupture  nucleated  in  a  localized 
region  of  high  stress,  the  average  stress  drop  over  the  rupture  area  might  be 
much  lower  than  the  initial  stress  drop,  while  radiated  energy  would  be  low  due 
to  the  gradual  stopping.  Considering  the  relative  uncertainties  of  these  cal¬ 
culations,  these  results  are  in  reasonable  agreement. 


TABLE  1 


Spectral  Parameters 


0153 

0356 

•13  +  .04  cm. sec 

•25  +  .07  cm* sec 

2.2  +  .6  cm2/sec 

3.1  +  .8  cm^/ sec 

1.2  h  z 

•  8  hz 

5.1  hz 

3.6  hz 

•24  sec 

•37  sec 

TABLE  2 


Rupture  DuraHnn  £ 
v 


Method 

0153 

0356 

Corner  Frequency 

•41  +  .1  secs 

•  62  +  .2 

Characteristic  Frequency 

•37  +  .1  secs 

•53  +  .15 

Rise  Time 

•39  +  .08  secs 

•60  +  .1 

Average 

•39  +  .05  secs 

•59  +  .1 

secs 

secs 

secs 

secs 


TABLE  3 


Source  Parameter 


0153  Event 

Radius  -  a  ■  1.2  ka 

Moment  -  m  •  3.5  +  .8  x  102**  dyne* cm 
o  — 

Stress  Drop  -  Ao  «  890  bars  Range  -  600-1100  bars 

Effective  Stress  -  x  ■  1040  +  350  bars 

e  — 

Radiated  Energy  -  Ea  -  8.7  +3. Ox  1020  dyne* cm 

Apparent  Stress  -  x  160  +  60  bars 
a  — 


0356  Event 

Radius  -  a  ■  1.65  km 

Moment  —  m^  ■  6.7  +  1.5  x  102lf  dyne* cm 

Stress  Drop  -  Ao  -  650  bars  Range  -  350-800  bars 

Effective  Stress  -  x  -  780  +  250  bars 

e  *— 

Radiated  Energy  —  E  ■  12.4  +  4.0  x  102®  dyne*cm 

Apparent  Stress  -  x  -  120  +  50  bars 
a  — 


100 


o-/pvv 

SGB 

SDP 

SQH 

NGI 

CNB 


cm/slc  _  (oi/sec)2 


cm/sec  (cm/sec)2 


log(displacement  amplitude) 


0356  Event 


MODEL  VARIATION  OF  SHEAR  WAVE  V 
CIRCULAR  MODELS 


PULSE  SHAPES  (e  =  30°) 

P 

ASYMMETRICAL  MODELS 
(v  =  .60) 


rupture  away  from  observer 


rupture  towards  observer 


BODY  WAVE  ANALYSIS  OF  THE  ST.  ELIAS  EARTHQUAKE 


In  order  to  Investigate  the  complexity  of  this  earthquake, 
it  is  useful  to  consider  both  P  and  S  body  waves  in  the  period 
range  from  2  secs  to  50  secs.  Using  three  WWSSN  long-period  seis¬ 
mograms  and  the  1-75  Benioff  seismograms  recorded  at  Palisades, 
we  have  established  that  the  earthquake  was  made  up  of  three 
distinct  sub-events,  preceded  by  a  small  initial  event.  The 
direction  of  rupture  propagation  appears  to  be  to  the  southeast. 

In  teleseismic  body  waves  from  a  shallow  fault,  waveform  com¬ 
plications  occur  through  the  Interference  of  the  depth  phases 
(i.e.,  the  phases  reflected  from  the  free  surface)  and  the  di¬ 
rect  phase.  To  bake  out  this  interference,  we  construct  a  free- 
surface  operator,  FSO(t),  for  each  body-wave  arrival,  using  the 

appropriate  radiation  patterns  and  reflection  coefficients  for 
T o iv  p^tTe-^wj  owe  'V. 

each  phase  to  determine  its  relative  amplitude. ^ Th  operator  is 

then  deconvolved  from  the  actual  pulse  shapes  in  order  to  obtain 

approximate  whole  space  (AWS)  pulse  shapes,  i.e.,  pulse  shapes 

without  the  interference  effects  of  the  free  surface.  The  delays 

of  the  depth  phases  were  adjusted  by  minimizing  the  deconvolutional 

noise,  thereby  obtaining  an  average  source  depth.  This  approach 

presumes  that  the  pulse  shapes  of  the  depth  phases  are  identical 

to  the  pulse  shapes  of  the  direct  phases.  Since  far-field  pulse 

shapes  depend  only  on  the  vector  slowness  of  the  body  wave  along 

the  fault  plane,  for  a  horizontal  fault  plane  the  pulse  shapes 


are  identical.  Thus  this  analysis  is  well  suited  to  earthquakes 
occurring  on  very  shallowly  dipping  fault  planes. 

The  steps  of  the  analysis  are  shown  in  Figure  1,  using  the 
F-wave  arrival  at  station  HKC  (Hong  Kong)  as  an  example.  The  lower 
most  trace  is  the  bandpassed  long-period  seismogram;  the  trace 
above  it  is  the  deconvolved  (from  the  instrument  response)  ground 
velocity.  For  this  arrival,  FSO(t)  has  a  positive  pulse  for  the 
direct  P,  a  very  small  positive  pulse  for  the  pP  and  a  large  nega¬ 
tive  pulse  for  the  s? •  The  convolutional  inverse,  labelled  IFSO(t) 
is  shown  next  to  it.  The  result  of  the  deconvolution  is  the  AWS 
velocity,  which  is  then  Integrated  to  obtain  the  AWS  displacement, 
shown  at  the  top  of  the  figure.  The  AWS  pulse  shapes  clearly  show 
the  three  major  sub-events,  which  are  labelled  1,  2,  and  3.  The 
initial  event  shows  up  as  a  small  step  in  the  seismogram  which 
becomes  a  single  bump  in  the  velocity^  tvcitA  - 

The  AWS  pulse  shapes  for  the  five  body  wave  arrivals  are  com¬ 
piled  in  Figure  2  along  with  their  epicenter  to  station  azimuths. 
There  is  substantial  variation  with  take-off  direction  and  wave- 
type.  The  pulse  shapes  at  azimuths  away  from  the  direction  of 
rupture  propagation,  i.e.,  the  HKC  P-wave  and  the  KEV  S-wave,  show 
the  longest  pulse  rise  times  (i.e.,  the  time  from  the  onset  of 
event  1  to  the  peak  of  event  3) ,  while  the  PAL  S-wave  has  the 
shortest.  Note  that  because  these  pulse  shapes  depend  on  the  slow¬ 
ness  of  the  body  wave,  the  S-wave  pulse  shapes  vary  more  strongly 
than  the  P-vave  pulse  shapes.  This  accounts  for  some  of  the  dlf- 
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ference  between  the  P-wave  and  the  S-wave  at  Palisades.  The  con¬ 
structive  interference  of  the  rupture  propagation,  as  seen  in  the 
S-wave  pulse  shape,  smooths  the  three  sub-events  into  a  single 
pulse.  The  marked  separation  of  events  2  and  3  in  the  PAL  P-wave 
may  be  the  result  of  a  slight  difference  in  the  focal  mechanism 
of  the  two  events,  as  this  arrival  is  at  a  P-pP-sP  node. 

Using  the  pulse  rise  times,  ,  of  the  sub-events,  we  can 
calculate  their  rupture  lengths  from  the  relation 

%  "  v  (1  "  Z  cos  c) 

where  l  is  the  rupture  length,  C  is  the  angle  between  the  direction 

,  of  v/eUaci 

of  rupture  propagation  and  the  takeoff  direction  of  the  body  wave 

and  v  is  the  rupture  velocity,  assumed  to  be  2.5  km/sec.  This 

calculation  gives  i  ■  9,  2A  and  16  km  for  the  three  events.  The 

PAL  S-vave  gives  a  total  rupture  length  of  68  km;  however,  as 

the  average  rupture  velocity  for  the  whole  event  must  necessarily 

be  smaller  than  the  rupture  velocities  of  the  sub-events,  this  is 

an  overestimate. 

The  average  rupture  depth  was  determined  using  the  PAL  pulse 

shapes,  as  the  free-surface  operators  for  these  arrivals  were  the 

most  sensitive  to  the  delay  times  of  the  depth  phases.  A  sS  delay 

of  6  seconds  and  a  sP  delay  of  5  seconds  gave  the  least  deconvolu- 

5«uru. 

tional  noise,  fixing  the  average  swptwre  depth  at  *  11  km. 
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The  deconvolutional  analysis  also  aids  Che  determination  of 
the  body-wave  moment,  as  it  coalesces  the  amplitude  information  of 
the  direct  phase  and  the  depth  phases  into  a  single  pulse.  The 
moments  of  the  sub-events  were  calculated  to  be  .8,  3.5  and 
7.6  *  1026  dyne  -cm,  respectively,  so  that  the  total  body-wave 
moment  is  1.2  x  1027  dyne-cm.  The  initial  event  has  a  moment 
<  4  x  1025  dyne  -cm.  Because  the  pulse  shapes  for  events  2  and  3 
are  not  separated  on  most  of  the  arrivals,  the  division  of  the  cumu¬ 
lative  moment  between  these  events  is  somewhat  ambiguous. 


John  Boatwright 

Lamont-Doherty  Geological  Observatory 
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FIGURE  CAPTIONS 


Figure  1. 


Figure  2, 


Figure  3. 
(optional) 


Steps  of  the  deconvolutional  analysis,  as  applied  to 
the  HKC  P-wave.  The  lowermost  trace  is  the  bandpassed 
(1.2  secs  to  66  secs)  long-period  seismogram.  The 
trace  above  it  is  the  deconvolved  ground  velocity, 
which  has  been  corrected  to  extract  the  effect  of  the 
high  pass  filter  to  a  unipolar  displacement  pulse. 

The  FSO(t)  and  IFSO(t)  operators  are  described  in  the 
text.  The  AWS  velocity  is  the  result  of  deconvolving 
FSO(t)  from  the  ground  velocity.  The  time  marks  are 
at  100  sec  intervals. 

Approximate  whole  space  displacement  pulse  shapes, 
obtained  from  the  body-wave  arrivals,  at  four  stations. 
$  is  the  azimuth  from  the  epicenter  to  the  station, 
as  measured  from  north.  The  numbers  identify  the 
sub-events  discussed  in  the  text.  The  dashed  line 
on  the  PAL  P-wave  is  the  presumed  baseline,  as  the 
sample  could  not  be  adequately  filtered  to  obtain  a 
flat  baseline. 

Radiation  patterns  for  the  phases  of  the  P  and  SH 

arrivals  analyzed.  The  squares  are  the  takeoff  angles 

TVsift  o-re 

to  the  stations;  from  left  to  right  ^  HKC,  KEV,  ESK, 
and  PAL.  The  sP  angles  have  been  corrected  for  the 


S-P  reflection  at  the  free  surface 


ELEMENTARY  SOLUTIONS  TO  LAMB'S  PROBLEM  FOR  A  POINT  SOURCE  AND  THEIR  RELEVANCE 


TO  THE  STUDY  OF  SPONTANEOUS  CRACK  PROPAGATION  ’U  THREE  DIMENSIONS 


By  Paul  G.  Richards ^ 


ABSTRACT 


Certain  exact  solutions  to  Lamb's  problem  (the  transient  response 
of  an  elastic  half-space  to  a  force  applied  at  a  point)  involve  the  com¬ 
putation  merely  of  three  square  roots,  and  about  ten  arithmetic  operations 
(+,  -,  x,  t) .  They  arise  when  both  source  and  receiver  lie  on  the  free 
surface.  It  is  just  these  solutions  which  are  needed  in  a  method  due  to 
Hamano  for  obtaining  the  slip  function  (displacement  discontinuity) ,  as 
a  function  of  space  and  time,  for  planar  tension  cracks  and  shear  cracks 
which  grow  spontaneously  with  arbitrary  shape.  The  solutions  are  described 
here  in  detail,  for  an  elastic  medium  with  general  Poisson's  ratio.  They 
include  perhaps  the  simplest-possible  example  of  the  P  wave. 


Lamont-Doherty  Geological  Observatory  Contribution  Number  0000. 

Lamont-Doherty  Geological  Observatory  and  Department  of  Geological 
Sciences  of  Columbia  University,  Palisades,  New  York  10964 


INTRODUCTION 


A  thorough  understanding  of  motions  in  an  elastic  half-space,  sub¬ 
jected  to  an  applied  force,  is  an  essential  part  of  wave-propagation 
theory  needed  to  interpret  seismic  waves.  For  this  reason,  half-space 
problems  have  been  the  subject  of  an  enormous  literature,  beginning  with 
Lamb's  (1904)  classic  study  of  displacements  set  up  by  forces  applied  at 
a  point  and  along  a  line  on  the  free  surface.  Here,  I  give  some  new 
solutions,  these  being  the  horizontal  motions  of  the  free  surface  for  a 
horizontal  force  applied  as  a  step  in  time  at  a  point  also  in  the  free 
surface.  (Throughout  this  paper,  the  half-space  is  oriented  with  a  hori¬ 
zontal  free  surface.  Taking  cartesian  axes  with  as  the  depth  coordinate, 
into  the  half-space,  the  free  surface  is  x3  *  0.) 

Whatever  the  value  of  Poisson's  ratio,  the  new  solutions  (which  aug¬ 
ment  the  work  of  Pekeris,  1955;  Chao,  1960;  and  Mooney,  1974)  are  ex¬ 
tremely  simple  to  compute.  However,  these  formulas  would  be  only  a  minor 
curiosity  if  it  were  not  for  one  very  important  application,  in  which 
speed  of  computation  is  essential.  This  application  is  suggested  by 
Hamano's  (1974)  method  for  studying  spontaneous  crack  propagation.  Since 
it  is  the  larger  problem  of  crack  propagation  which  has  motivated  the 
present  study,  I  shall  in  what  follows  give  a  brief  review  of  Hamano's 
method,  before  giving  the  simple  solutions  to  Lamb's  problem. 
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MOTIVATION 

Within  an  infinite  homogeneous  elastic  medium,  initially  at  rest, 
suppose  that  a  crack  nucleates  at  time  t  =  0  and  subsequently  grows 
within  the  plane  x3  =  0.  Then  a  useful  representation  of  displacement 
u  =  u(x,t)  throughout  the  medium  can  be  given  as 

OO  00  00 

un  (x,t)  -  /  dT  f  f  dCjdC2  G  (x,t-T;C  5  0,0)  T  (5  ,5  ,t) .  (1) 

—oo  —oo—oo  r  -  I  pit. 

Here,  G  (x,t;£,T)  is  the  Green  function  for  the  medium,  being  the  n- 
component  of  displacement  at  (x,t)  due  to  a  unit  impulse  applied  in  the 
p-direction  at  position  £  and  time  t.  For  purposes  of  computing  G^ , 
it  is  required  that  the  whole  plane  x^  =  0  be  a  traction-free  surface. 

T  (£j,£2,t)  is  the  actual  traction  occurring  on  the  whole  plane  x3  =  0  of 
which  the  crack  is  a  part.  Equation  (1)  is  described  further,  and  proved, 
by  Das  and  Aki  (1977a)  and  Aki  and  Richards  (1979,  their  equation  2.43).  In¬ 
tuitively,  the  above  representation  can  be  understood  as  replacing  the 
actual  (crack)  source  of  radiation  by  a  whole  plane,  separating  the  medium 
into  two  half-spaces.  Into  each  half-space,  waves  are  radiated  due  to 
the  same  tractions  as  those  set  up  by  the  crack,  applied  over  the  half¬ 
space  surfaces.  From  the  space-time  element  dxd£  d£  there  is  an  applied 
impulse  of  strength  dtd^d^  T^  (£j,£2,x)  in  the  p-direction.  If  the 
displacement  contribution  to  un  (x,t)  from  this  element  is  to  be  considered 
in  isolation  from  tractions  acting  elsewhere  on  X3  =  0,  then  the  ap¬ 
propriate  Green  function  Gnp  must  be  constrained  by  having  zero  traction 
over  the  surface  of  the  half-space.  Hence,  it  must  be  a  solution  to  Lamb’s 
problem. 
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Hamano  (1974)  pointed  out  that  for  shear  cracks  and  tension  cracks, 
a  soluble  scheme  for  the  displacement  discontinuity  [u],  say,  across  the 
crack  can  be  set  up  from  equation  (1)  by  considering  the  x  position  it¬ 
self  in  the  crack  plane,  x  =  (xj,x2,0),  and  using  symmetry  properties  of 
u  and  T  across  x3  =  0  to  constrain  the  displacement  and  traction  on  dif¬ 
ferent  parts  of  the  plane  of  the  crack. 

For  a  general  tension  crack,  [u]  =  [0,0, U3]  and  T  =  (0,0, T3)  so  that 

the  only  Green  function  needed  is  G33. 

For  a  general  shear  crack,  [u]  =  [uj,u2,0]  and  T  =  (Tj,T2,0)  so  that 
the  only  Green  functions  needed  are  G^,  G12,  G21  and  G22  (see  Das  and  Aki, 
1977a,  for  a  related  study  of  two-dimensional  cracks).  The  jump  in  u3  is  zero, 
because  opposite  faces  of  the  crack  remain  in  contact.  Tj  is  zero,  because 
planar  shearing  cannot  change  the  normal  stress  on  the  crack  plane.  Together 
with  G33  for  tension  cracks,  these  five  different  Lamb  problems/Green  functions 
need  be  studied  only  for  the  case  that  both  source  and  receiver  lie  in 
the  free  surface  of  the  half-space.  Hamano' s  method  is  important  in 
offering  the  chance  to  study  spontaneous  crack  growth  for  completely  gen¬ 
eral  si  apes  of  planar  cracking.  This  paper  contributes  to  that  goal,  by 
showing  that  just  these  Green  functions  are  almost  trivially  simple  to 
compute.  For  completeness,  a  single  integral  is  also  given  below,  in  terms 
of  which  the  remaining  four  components,  G13,  G23,  G31,  G32,  can  efficiently 
be  computed. 


FORMAL  STATEMENT  OF  PROBLEM,  AND  ITS  SOLUTION 
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In  this  section,  explicit  formulas  are  derived  for  Gnp  (x1,x2,0,t;C1,C2,0,0) , 
this  being  the  n-  component  of  displacement  at  position  (xj,x2,0)  and  time  t, 
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within  the  free  surface  (x3  =  0)  of  a  homogeneous,  isotropic  elastic  half¬ 
space,  due  to  a  unit-step  force  in  the  p-direction  applied  also  within 
the  free  surface,  at  position  (Sj,^*®)*  the  step  occurring  at  time  0. 

Once  this  solution  has  been  found,  G  for  an  impulse  (as  required  in  the 
section  above)  is  given  by 


(x,t-t;C,0) 


3_  G  H 
3  s  np 


(x,s;£,0) 


s 


t-T 


Since  related  problems  have  had  such  a  wide  exposure,  I  shall  abbre¬ 
viate  the  description  of  how  a  solution  is  obtained.  Thus,  in  general, 
the  solution  to  Lamb's  problem  for  a  point  source  can  be  obtained  as  an 
integral  over  just  one  variable.  In  our  case. 


H  ^  ^  HfT-13  L  P  dP 

Gn  (x  x,,0,t;0,0)  =  *2  ImagU  -V  -  } 

1  (T2-P2)  2  [ (A-2P2) 2+4XYP2 ] 

where  p  =  rigidity,  r  =  (x12+x22)  2»  H  is  the  unit  Heaviside  step  function, 
and  capital  letters  in  the  integrand  denote  dimensionless  quantities: 

A  =  ct2/62(a  =  P-wave  speed,  8=  S-wave  speed) 

T  *  at/r  (T  =  1  being  the  P-wave  arrival  time) 

X  ■=  (l-P2)^  or  -i  (P2-!)*5  ,  Y  =  (A-P2) ^  or  -i  (P^A)*5, 

Imag  {}  denotes  the  imaginary  part  of  {},  and 
Ln  =  {(T2-P2)  [2Y-4X+(A-2P2)/Y]  +  AY}  COS2<() 

-  Cr2[2Y-4X+(A-2P2)/Y]  -  AY)  sin2<f> 

Li2  ■  L21  ■  (2T2-P2)  [2Y-4X+(A-2P2)/Y]  cos*  sin* 


(1) 
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1*22  -  {(T2-P2)[2Y-4X+(A-2P2)/Y]  +  AY)  sin2* 

-{ T2 [ 2Y-4X+(A-2P2) /Y ]  -  AY)  cos2* 

and 

L13/COS*  »  1,23/sin*  -  -  L31/cos*  *  -L32/sin*  *  T(A-2P2)  -  2TXY,  where  * 
is  given  by  x1  »  rcos*,  x2  =  rsin*,  so  that  *  is  the  azimuth  to  x. 

Formula  (1)  can  be  written  down  from  Johnson  (1974,  his  equations  26- 
34,  but  using  P2  for  his  a2t2r-2  -  a2p2).  Our  variable  P  is  a  (dimension¬ 
less)  horizontal  slowness,  and  (1)  is  essentially  a  Cagniard  solution  in 
the  form  advocated  by  Helmberger  (1968) .  Both  source  and  receiver  lie  in 
the  free  surface,  so  the  Cagniard  path  lies  just  above  the  real  P-axis  in 

the  complex  P-plane  (see  Figure  1,  and  caption). 

H  H  H  H  II 

In  fact,  the  P-integrals  for  Gn  ,  Gi2  »  G2i  ,  G22  and  G33  can  be 
given  in  closed  form.  This  is  not  possible  for  the  four  remaining  entries 

17 

in  G  ,  but  these  four  are  all  proportional  to  just  one  integral  which  is 
still  fairly  simple  to  compute.  We  note  first  that 

GnH  (x!  ,x2,0,t;0,0)  =  [Ij  (T)  cos2*  -  I2(T)  sin2*]/(irpr) 
g12H  *  G2lH  *  [I1+I2]  cos*  sin*/(irpr) 

G22H  *  [Xi  sin2*  -  I2  cos2*]/(iryr)  ^ 

G33H  -  I3(T)/Oryr) 

where  arguments  have  been  written  out  explicitly  only  for  the  first  of 
equations  (3) .  Just  three  functions  of  dimensionless  time  are  needed  to 

17 

evaluate  these  five  components  of  G  ,  namely 
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Il(T) 


i2(t) 


iImag  {/  JL(IZ1^,  .CI.2-P2)[2Y-4X+(A-2P2)/Y]  +  AY  pdp  } 
1  (T2-P2r  (A-2P2) 2+4XYP2 


—  Imag  {  / 


T  H(T-l)  42[2Y-4X+(A-2P2)/Y1  -  AY 


(T2-P2> 


(A-2P2)2+4XYP2 


PdP  }  (4) 


and 


I3(T)  -  J  Ima?  {  /  -(— 1)-u  - — - 

1  (T2-P2)  2  (A-2P2) 2+4XYP2 


PdP  }  . 


It  follows  from  (3)  that  Ij  and  I3  are  displacements  within  the  vertical 
plane  containing  source  and  receiver,  like  the  dominant  motion  in  P-SV,  whereas  I2 
is  a  displacement  transverse  to  this  plane,  like  the  dominant  motion  in  SH.  I 
give  analytic  expressions  below  for  each  of  these  three  integrals.  A  fourth 
dimensionless  solution  to  Lamb's  problem  is  introduced  via 


Gl3H  (xi,x2,0,t;0,0)  *  Ii+(T)  cos<(>/ (iryr)  , 

U  H  u 

G23  “  I4  sin<J>/(iryr) ,  G31  =  -  I4  cos<t>/(Tryr) ,  and  G32 

The  integral  for  this  solution  is 

I„(T>  -  i  pdp  ) 

1  (T2-P2K  v  (A-2P2)2  +  4XYP2, 


(5) 

-  I4  sin<f>/(irijr)  • 


(6) 


which  cannot  be  given  in  closed  form. 
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For  a  vertical  point  force,  Pekeris  (1955)  obtained  a  closed-form  solu¬ 
tion  for  the  vertical  displacement  (our  I3)  and  a  sum  of  elliptic  integrals 
for  the  horizontal  displacement  (our  I4) .  His  solutions  were  restricted 
to  the  case  a2/$2  =  3  (Poisson's  ratio  =  0.25),  and  Mooney  (1974)  indicated  how 
the  evaluation  of  I3  might  be  carried  out  for  any  value  of  a2/(32  (though  Mooney 
did  not  publish  the  solution  formulas) .  For  a  horizontal  point  force,  in  a 
medium  with  a2/02  ■  3,  Chao  (1960)  obtained  closed-form  solutions  for  the  hori¬ 
zontal  displacements  (our  1^  and  I2)  and  a  sum  of  elliptic  integrals  for  the 
vertical  displacement.  Solutions  themselves  have  not  previously  been  given 
explicitly,  for  general  a2/B2,  for  any  one  of  the  four  basic  (dimensionless) 
solutions  1^. 

In  every  case,  the  basic  approach  involves  writing 


(A-2P2) 2+4XYP2  (A-2P2)4  -16X2Y2P4  Cubic (P2) 

so  that  the  new  denominator,  of  sixth  order  in  P,  is  real  throughout  the 
P-axis  integration  and  has  no  branch  cuts.  The  imaginary  parts  of  the  new 
Integrals  are  easy  to  identify  (together  with  a  semi-residue  contribution  to 
I4  from  indenting  around  the  Rayleigh  pole:  see  Figure  1  caption).  If  roots 
Rj,  R2 ,  R3  are  found  for  the  cubic  in  P2,  integration  for  I*,  I2,  I3  becomes 
possible  using  the  partial-fraction  decomposition 

- 1 -  ,_JL.  +  _b_  +  _£_ 

(A-2P2) 4-16X2Y2P2  P2-Rj  P2-R2  P2-R3 


(8) 
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The  Rayleigh  pole  lies  at  P  *  R32  *  a/y  (always  on  the  real  P-axis, 
just  to  the  right  of  P  =  0/6,  since  Rayleigh  wave  speed  y  is  a  few  percent 
less  than  g) .  If  Poisson's  ratio  is  less  than  a  critical  value,  approx¬ 
imately  0.263,  then  Rj  and  R2  are  real  and  lie  between  0  and  1.  But,  for 
greater  values  of  Poisson's  ratio,  Rj  and  R2  are  complex  conjugates,  as 
are  a  and  b  in  (8),  though  c  is  always  real.  In  this  case  there  are  poles  in 
the  complex  P-plane  which  can  be  associated  with  the  so-called  F  wave  (Gilbert 
st  al.t  1962;  Chapman,  1972;  Aki  and  Richards,  1979),  appearing  between  P- 
and  S-arrivals.  Rj  and  R ^  lie  on  a  Riemann  sheet  different  from  that  which 
contains  the  Cagniard  path. 

Substitution  of  (8)  into  (7)  and  (4)  leads  to  24  integrable  terms  for 
each  of  I j  and  1 2*  and  6  such  terms  for  1 3.  Extensive  cancellation  does 
eventually  occur.  The  solutions,  involving  real  positive  constants 
Cj  (j  =1,  ...  ,  7)  and  complex  constants  (k  *  1,  2,  3),  are  as 
follows: 

For  times  prior  to  and  including  the  P-arrival,  T  _<  1, 

II  *  I2  "  13  "  (9) 

For  times  between  P-  and  S-arrivals,  1  <  T  <  a/ 6,  there  are  two 
different  kinds  of  elastic  media  to  consider.  If  Poisson's  ratio  is 
less  than  0.263  (corresponding  to  A  -  a2/g2  <  3.11), 

11  =  T2[Cl  (T2^)**5-  c2(T2-R2)~^  -  c3 (R3-T2)”**] 

12  -  -c4  +  Cl  (TZ-Rj)*  -  c2  (T2-R2),s  +  c3  (Rj-T2)1* 

1 3  ■  c4  -  c5  (T2-Rx  )“Ss  +  eg  (T2-R2)->S  -  c7  (R3-T2)_>S  . 


(10) 
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If  Poisson's  ratio  is  greater  than  0.263,  it  is  necessary  first  to  define 
the  complex  square  root  (complex,  because  R2  is  then  complex), 

CROOT  -  [(1-R2)(T2-R2)]^  (11) 

in  which  the  choice  of  sign  is  made  such  that  the  complex  number 
1  +  2  (1-R2  -  CROOT) /(T2-l)  has  magnitude  less  than  unity.  Then, 

11  -  -T2  [Real  {C i / CROOT }  +  c3  (R3-T2)_ls] 

12  -  -C4  -  Real  (C2  *  CROOT)  +  c3  (R^T2)*5  (12) 

13  -  c4  +  Real  {C3 /CROOT)  -  C7  (R3-T2)-,s  . 

For  times  between  the  S-arrival  and  the  Rayleigh  wave  arrival, 
a/6  <  T  <  a/y, 

Ij  -  0.5  -  2c3  T2 (R3-T2)-^ 

12  -  -2c4  +  2c3  (Ra-T2)*5  (13) 

13  -  2c4  -  2c7  (R3-T2)_Js 

For  times  after  the  Rayleigh  arrival,  a/y  <  T, 
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11  -  0.5 

12  -  -2C4  (14) 


I3  *  2c4  • 


Constants  In  the  above  solution  need  be  evaluated  just  once  for  a 
given  elastic  medium,  specified  by  the  ratio  a2/S2.  An  effective  approach 
Is  first  to  find  the  largest  root  R3  of  the  Rayleigh  cubic;  then  to  factorise 
P2-R3  from  the  cubic  and  solve  a  quadratic  for  Ri  and  R2.  Constants  a,  b,  c 
in  (8)  are  given  by 


a"1  -  16(A-l)(R1-a2)(R3-R1) 


b"1  -  16 (A-l) (Rx-R2) (R2-  R3) 


(15) 


c'1  -  16(A-1)(R3-R1)(R2-R3)  . 


Then 


cj  *  -2aA(A-Rx) (l-Rx)*4 


c2  -  2bA(A-R2)(l-R2)Js 


c 3  -  -2cA(R3-A)(R3-l)Js 


C4  »  A/(8A-8) 


c5  — 2aARx(l-Ri)(A-Rx) 


% 


c6 


2bAR2 (1-R2) (A-R2) 


H 
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u 

c7  -  -2cAR3(R3-1)(R3-A)  d  =  4bA(A-R2) (1-R2) 

C2  -  4bA(A-R2)  C3  =  bA(A-2R2)2(l-R2). 

Solutions  given  above  for  1^,  I2,  I3  require  at  most  the  evaluation  either 
of  three  real  square  roots,  or  (depending  on  Poisson's  ratio)  the  eval¬ 
uation  of  one  complex  square  root  and  one  real  square  root.  These  (worst) 
cases  occur  only  between  P-  and  S-arrivals.  In  terms  of  these  closed-form 
solutions,  all  the  five  components  of  G  relevant  to  Hamano's  method  for 
studying  spontaneous  shear  and  tension  cracks  can  be  rapidly  computed  via  (3) . 

Although  I4  can  be  given  in  terms  of  elliptic  integrals  (with  complex 
arguments  when  Poisson's  ratio  is  greater  than  0.263),  it  is  probably  more 
efficient  directly  to  integrate  as  follows: 


0  for  T  <  1; 

—  /  ^~AH4l.p2) ^(A ~2p2-I  dx  for  1  <  T  <  a/8, 

*  0  (A-2P2)4-16X2Y2P4 


I4(T)  -1)  where  P2  *  (T2-l)  sin2X  +  1; 

i 

2TA  J/Z  (P2-l) (A-P2) (A-2P2)dV 
17  0  (T2-P2)  (A-2P2)  4-16X2Y2P4  ] 

\ 

where  P2  *  (A-l)sin2Y  +  1. 

\ 


H<T-^>V  tor  a/8  <  T, 

(T2-R  )** 

3 


(16) 


Integrals  with  respect  to  X  and  V  here  have  well-behaved  integrands.  Note 
that,  at  time  T  -  a/y  ■  R3'2,  a  singular  Rayleigh  wave  arrives  (see  Figure  lb 
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caption)  with  strength  proportional  to  the  positive  real  constant 

c  ■  *5cA(A-2R3)  3/R3 .  (17) 

8 

In  Figure  2  are  shown  the  time-dependences  of  Ii,  I2,  I3  for  four 
different  values  of  a2/ 82.  We  note  the  following  basic  properties:  (i) 
Displacements  Ij  and  I3  are  continuous  across  the  P-arrival,  as  are  1 2  and 
the  particle  velocity  dl2/dT.  These  results  follow  from  (10),  (12),  and  rela¬ 
tions  between  constants  appearing  in  these  formulas,  (ii)  Ij  and 
I3  are  continuous  across  the  S-arrival,  but  have  discontinuous  slopes, 
whereas  I2  itself  is  discontinuous,  (iii)  I 2  is  continuous  across  the 
Rayleigh-wave  arrival  time,  but  1^  and  I3  are  singular.  All  three  solu¬ 
tions  are  exactly  constant  after  the  Rayleigh  singularity:  these  constants 
must  then  be  the  static  solutions,  (iv)  For  the  horizontal  displacement  due 
to  a  horizontal  force,  the  step  (in  I2)  at  the  S-arrival  can  be  seen  from 
(3)  as  having  the  orientation  of  an  SH-wave,  whereas  the  singularity 
(in  Ij)  at  the  Rayleigh-arrlval  occurs  as  P-SV  motion.  However,  because 
P-wave  motion  is  not  in  general  exactly  longitudinal,  the  transverse  motion 
(given  by  I2)  does  begin  at  the  P-wave  arrival,  (v)  A  P  wave  is  apparent 
in  I3  at  times  between  and  P-  and  S-arrivals,  becoming  more  apparent  with 
Increasing  values  of  o2/62.  Since  it  arises  from  a  single  algebraic 
expression,  the  term  Real  {C3/CROOT}  in  (12),  detailed  properties  of  this 
wave  are  easy  to  investigate. 

In  Figure  3  is  shown  the  time-dependence  of  I4  for  four  different 
values  of  a2/S2.  Romberg  integration  was  used,  requiring  occasionally  up 
to  128  intervals  for  1%  accuracy.  There  is  a  discontinuous  slope  at  the  P-  and 
S-arrivals;  a  jump  to  a  singularity  at  the  Rayleigh  arrival;  and  thereafter 
a  gradual  decay  to  the  static  limit. 
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CONCLUSIONS 


Perhaps  the  main  achievement  of  this  paper  Is  the  exact  form  of  constants 

c  ,...., c  (positive  real.  If  used)  and  complex  constants  C  ,  C  ,  C  ,  in  terms 
18  12  3 

of  which  the  complete  solution  to  Lamb's  problem  can  be  given  for  any  orient¬ 
ation  of  applied  force,  any  displacement  component  and  any  value  of  Poisson's 
ratio,  provided  both  source  and  receiver  lie  in  the  free  surface. 

Four  scalar  solutions  in  G,  involving  the  cross-terms  (vertical  or  hori¬ 
zontal  displacements  due  respectively  to  horizontal  or  vertical  applied  force) , 
cannot  be  given  in  closed  form,  but  a  well-behaved  integral  solution  is 
possible  in  general. 

In  the  case  of  horizontal  displacements  due  to  a  horizontally  applied 

force,  the  solutions  are  relevant  to  a  method  for  studying  spontaneous  shear 

cracks.  For  the  case  of  vertical  displacement  due  to  a  vertical  force,  the 

solution  has  relevance  to  tension  cracks.  In  both  these  cases,  solutions 

(3),  (9)-(14),  for  a  step-applied  force,  are  so  simple  that  the  following 

can  readily  be  derived  in  closed  form:  (a)  solutions  for  an  impulsively- 

applied  force;  (b)  solutions  averaged  over  (r,  r+Ar),  (<(> , <|>+A<J> )  and  (t,t+At); 

(c)  14  of  the  displacement  fields  3G  /3£  *  G  due  to  a  single-couple. 

np  q  np , q 

Specifically,  we  can  use  reciprocity  on  G  so  that  the  derivative  is  con¬ 
ducted  with  respect  to  receiver  coordinates.  From  the  five  closed-form 
solutions  in  G,  10  single-couple  displacement  fields  can  be  obtained  by 
differentiating  in  the  1-  and  2-directions,  parallel  to  the  free  surface. 

The  solutions  for  (n,p,q)  ■  (1,3,3),  (2,3,3),  (3,1,3)  and  (3,2,3)  can  also 
be  recovered  in  closed  form,  by  using  the  linear  strain  constraints  at  a 


stress-free  surface. 
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The  simplicity  of  the  five  scalar  solutions  in  G,  which  are  associated 
with  Hamano's  method  for  studying  cracks,  is  so  remarkable  that  it  gives 
high  hopes  of  successful  development  of  a  3-dimensional  study  of  spontaneous 
shear  and  tension  fractures. 
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FIGURE  CAPTIONS 


The  Cagniard  integration  path  for  (1) ,  (4)  and  (6)  is  shown  as  a 
solid  heavy  line.  Singularities  of  these  integrands  are  shown 
in  (a).  They  consist  of  branch  cuts  trending  to  the  right  along 
the  real  P-axis  from  points  1  and  a/6,  and  a  pole  at  a/y  (y  being 

the  Rayleigh  wave  speed).  As  well  as  the  pole  at  a/y  »  R3 2,  we 

i-  V 

also  show  schematically  the  poles  at  R  2,  R2  2  in  the  right  half¬ 
plane.  Ri,  R2  and  R3  are  roots  of  the  Rayleigh  cubic  in  P2, 
(A-2P2)4-16X2Y2P4  =  16(1-A)(P2-R1)(P2-R2)(P2-R3).  For  A  <  3.11, 

R^  and  R2  are  real  and  lie  between  0  and  1.  But  for  A  >  3.11, 

Rj  and  R2  are  complex  conjugates.  Accurate  locations  for  different 
A  are  given  in  Figure  3b. 

In  (b)  is  shown  a  path  of  integration  in  the  vicinity  of  the 
Rayleigh  pole.  In  the  limit,  as  the  semi-circle  radius  shrinks 
to  zero,  integration  reduces  to  a  principal  value  integration 
plus  -in  x  residue  at  the  Rayleigh  pole.  The  residue  is  imaginary 
from  the  integrands  (4)  of  Ij,  I2,  I3,  and  hence  gives  zero  net 
effect  when  the  imaginary  part  is  evaluated  after  multiplication 
by  -in.  But  the  residue  from  integrand  (6)  for  I4  is  real, 
leading  to  a  non-zero  contribution  from  the  Rayleigh  pole  when 
T  >  a/y. 

Here  are  shown  the  fundamental  time-dependences  of  displacement 
for  (a)  Ij,  the  longitudinal  horizontal  displacement  for  a 
horizontally-applied  force;  (b)  I2,  the  transverse  horizontal 
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displacement  for  a  horizontally-applied  force;  and  (c)  I3,  the 
vertical  displacement  for  a  vertically-applied  force.  We  have 
chosen  to  plot  values  of  I3  positive  downwards,  so  that  upward 
values  in  (c)  correspond  to  -  I3,  and  hence  to  the  convention 
common  in  seismology  of  recording  vertical  motions  as  positive 
upwards. 

In  each  case,  the  time-dependence  is  worked  out  for  four 
different  values  of  a2/B2.  Time  for  these  four  cases  is  scaled 
so  that  P-arrivals  (at  T  »  1)  and  S-arrivals  (at  T  ■  a/B)  are 
aligned.  Values  are  plotted,  as  heavy  solid  lines,  only  between 
amplitudes  +1.  Ij  and  I3  in  fact  are  singular  at  the  Rayleigh 
arrival  (marked  R),  thereafter  jumping  immediately  to  the  static 
value. 

Dotted  lines  give  values  of  Ij,  I2,  I3  scaled  up  by  a  factor 
of  15,  and  hence  display  the  detailed  time-dependence  at  low 
amplitudes. 

Figure  3:  (a)  Values  of  the  fundamental  solution  I4  as  a  function  of 

time.  A  closed-form  solution  is  not  possible  in  this  case. 
Computation  is  for  four  different  ratios  of  ci2/B2.  Dotted  lines 
show  values  of  15  *  I4. 

(b)  Since  T  can  be  regarded  as  a  value  of  P,  the  dimensionless 

horizontal  slowness,  we  have  here  shown  the  complex  P-plane 

with  the  same  scale  as  the  T  axis  in  (a),  and  Figure  2a,b,c. 

u 

Values  of  15  *  I4  are  repeated  from  (a) .  Singularities  Rj  , 

R2  »  R3  here,  for  different  values  of  a2/B2,  occur  then  at  times 
(P-values)  which  are  indicative  of  what  turn  out  to  be  properties 
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of  the  P  and  Rayleigh  waves.  Thus,  R2  2  for  a2  *  2%(}2  is  almost 
coincident  with  the  ordinary  P-wave  arrival,  making  the  latter 
highly  impulsive  for  Ij  and  I3  because  of  the  term  in 

(T2-R2)~^.  At  larger  a 2/B2,  the  occurrence  of  complex  R2  2 
with  real  values  greater  than  one  leads  to  the  emergent  broad 
pulse  between  P-  and  S-arrivals  in  I3  and  I4.  It  is  interesting 
that  such  a  P  wave  is  not  apparent  for  Ij  and  I2. 
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